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NOMENCLATURE 


*  * 

(a  /ar,  vp/p  ),  speed  of  sound 

(also,  constant  in  burning  rate  law  r  =  apn  ) 

.  *.  *  . 

(A  /Ack'>  cross  sectional  area 

chamber  reference  area 

control  volume  perimeter  surface  area 

(5-space)  location  of  left-hand  boundary  of  flow  field 

region 

(B  /ay  ),  frequency  factor  in  Arrhenius  surface  reaction 
r  v  c 

(5-space)  location  of  right-hand  boundary  of  flow  field 
region 

specific  heat  at  constant  pressure  of  gaseous  combustion 
products 

specific  heat  of  solid  propellant  material 

specific  heat  at  constant  volume  of  gaseous  combustion 

products 

control  surface 

control  volume 

control  surface  perimeter 

discriminant  defined  in  Equation  (2.50) 

-ft  ft 

(e  /c  Tr),  specific  internal  energy  of  gas 
activation  energy  of  Arrhenius  flame  reaction 
activation  energy  of  Arrhenius  surface  reaction 
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h 
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hcon 
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Surface 
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* 
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net  reaction  force  acting  on  the  fluid  within  the 
control  volume 

represents  boundary  condition  Equation  (3 .20) (c )  at  time 
"n" 


(h  /CpTr>  p/p)  specific  enthalpy  of  gas 

enthalpy  of  combustion  products  entering  control  volume 

from  flame  zone 

convective  heat  transfer  coefficient 

thermal  conductivity  of  combustion  products 

thermal  conductivity  of  solid  propellant  material 
*  * 

(  2  )  »  constant  which  follows  from  Stokes  Flow 

PmCT  ar 


Drag  Law 

axial  reference  length 
*  * 

(m  / Prar ) j  propellant  surface  mass  flow  rate  per  unit 

surface  area 
.# ,  *  * 

(m  /p  a  mass  flow  rate  of  gas  per  unit  surface  area 
g'  r  r 
#  * 

(mp/prar),  mass  flow  rate  ci*  particles  per  unit  sui'face 
area 

(u/a),  Mach  number 
reaction  index 
(P  /Pr),  pressure 

heat  transfer  rate  tc  propellant  surface 
defined  in  Equation  (3.11) 
exothermic  heat  release  in  flame  reaction 
endothermic  heat  release  in  surface  reaction 
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sur 
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U,  U 

J  v 


(r  /aryc ),. propellant  surface  regression  rate 


value  of  r  when  no  disturbances  are  present 

*  i 

w 


*  ‘X* 

(R  /R  .  ),  local  chamber  radius 


(  -■  ^  ,  combustion  response  function 

V/5' 

the  real  part  is  computed  from  the  numerical  solution 

/  Am/m  \ 

Ap/p 

.  *  .  *. 

r.(BChA  ) 

R^  defined  in  Equations  (2.9)  through  (2.13) 

VA*h/n  ,  radial  reference  length 

-ft 

(Ro/W  ),  gas  constant 
universal  gas  constant 

entropy  of  gas,  defined  in  Equation  (2.l4) 

*  *  * 

(t  a^/L  ),  time 
(t/period  of  wave) 

'X*  *X- 

(T  /t^),  temperature 
defined  before  Equation  (3-23) 
adiabatic  flame  temperature 
flame  temperature 

temperature  of  cold  solid  propellant 
temperature  of  surface  of  burning  propellant 

-jf  -X* 

(u  /a  ),  velocity  of  gas 

(Up/a^),  velocity  of  solid  particles 

transformed  velocities  defined  in  Equations  (2.17)  and 

(2.18),  respectively 
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T',  see  Equation  (3-27) 

velocity  of  gas  normal  to  propellant  surface  in  flame 
shock  wave  velocity  relative  to  fixed  coordinates 
molecular  weight 

ft 

(x  /L  ),  axial  coordinate 

transformed  axial  coordinate  defined  in  Equation  (2.15) 

ft  ftA 

(y  /L  yc),  coordinate  normal  to  regressing  propellant 
surface 

.  #.  *  *  *  *  i 
<VVsarL. > 

location  of  edge  of  flame 

(«ATr> 

[(kgQ^)/(ksCpTrpsa;)]  »  (t) 

(m  cpL  yc/kg) 

defined  in  Equation  (2.l6) 
isentropic  index 
(d§/dx)/(c-b) 

defined  by  Equation  (3.32) 
thermal  conductivity  of  solid  propellant 
complex  function  of  frequency,  0 
(solution  of,  \(\  -  l)  =  in) 

5  =  ?(x),  stretched  axial  coordinate 


^1’  ^d  defined  Equation  (3-"^0) 


ft  ft 

(P  /pr)>  density  of  gas 

ft  ft 

(Pp/Pr )’  density  s°H-d  particles  in  chamber  volume 
density  of  solid  particle  material 


solid  propellant  density 
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Superscripts 
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Subscripts 
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average  particle  diameter 
(«)*L*/a*),  frequency 
reaction  rate  in  gas  phase  flame 
(to/r8) 

dimensional  quantity 
time-varying  fluctuation 

gas  phase 
solid  particles 
reference  state 
solid  propellant  material 


CHAPTER  I 


INTRODUCTION 

General  Introduction  and  Perspective 
The  full-thrust  operating  characteristics  of  solid  propellant 
rocket  motors  are  never  completely  time-independent.  The  performance 
of  the  system  is  acceptable  with  a  certain  amount  of  rough  combustion, 
defined  as  small  amplitude  unorganized  flow  field  oscillations.  On 
occasion,  larger  amplitude  organized  disturbances  appear  in  the  com¬ 
bustion  chamber  and  nozzle  flow  field  during  the  ignition  transient 
or  after  full-thrust,  has  been  obtained.  If  such  an  instability  reaches 
a  large  amplitude  limit  cycle,  it  may  lead  to  an  increase  in  mean  chamber 
pressure  and  burning  rate,  excessive  heat  transfer  rates,  and  a  severe 
vibration  level.  The  occurrence  of  any  one  of  these  may  result  in 
malfunction  or  destruction  of  the  engine. 

Organized  oscillations  are  evidence  of  a  positive  feedback  between 
the  propellant  combustion  process  and  the  flow  field  disturbance.  As 
noted  by  Price  ^ ,  a  rocket  motor  is  like  a  self-excited  oscillator 
where  combustion  is  the  energy  source,  the  flow  in  the  combustion 
chamber  and  nozzle  is  the  oscillator,  and  the  propellant  burning  rate 
sensitivity  to  flow  field  disturbances  is  the  coupling  between  the  two. 

A  combustion  instability  analysis  seeks  to  describe  how  this  interaction 
occurs,  determine  what  combination  of  system  parameters  may  lead  to  an 
unstable  situation,  and  suggest  ways  it  might  be  counteracted. 


Preceding  page  blank 
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In  general,  rocket  motor  instabilities  can  occur  in  longitudinal, 

radial,  or  transverse  modes  at  frequencies  which  are  close  to  the  natural 

acoustic  frequencies  of  the  combustion  chamber.  Many  solid  propellant 

rocket  motors  are  prone  to  undersirable  axial -mode  or  longitudinal 

(2) 

instabilities  in  an  intermediate  frequency  range  of  100  to  1000  Hz. 
The  infrequent  occurrence  of  high  frequency  modes  can  be  attributed  to 
the  use  of  metalized  propellants  which  produce  solid  particles  in  the 
gas  flow.  In  these  cases,  the  resultant  drag  between  gas  and  particles 
and  its  associated  energy  dissipation  contribute  a  large  damping  effect 
'to  the  unsteady  motion  in  the  chamber.  However,  experimental  obser¬ 
vation  and  analytical  considerations  indicate  that  particulate  damping 
is  ineffective  in  controlling  axial  mode  intermediate  frequency  insta¬ 
bilities.  In  certain  instances,  solid  rocket  engines  have  shown  an 

(3) 

increased  susceptibility  w/  to  axial -mode  instabilities  as  the  metal 
content  of  the  propellant  is  increased. 

Longitudinal  instabilities  normally  arise  spontaneously,  pre¬ 
sumably  triggered  by  the  presence  of  small  amplitude  unorganized  distur¬ 
bances  in  the  combustion  chamber,  or  by  a  momentary  flow  blockage  due 
to  a  solid  fragment  exhausting  through  the  nozzle.  Attempts  to  isolate 
controlling  parameters  and  establish  stability  trends  by  firing  full- 
scale  engines  become  prohibitively  expensive.  Outside  of  a  few  reported 
cases  of  instability  in  full-scale  engines,  most  experimental  data 

has  been  gathered  from  laboratory  ;ests  of  small-scale  motors.  A  very 

(7_12 ) 

extensive  investigation  has  been  carried  out  at  the  Canadian 

Armament  Research  and  Development  Establishment  (CARDE)  to  determine 

the  effect  of  chamber  diameter,  length,  throat-to-port  area  ratio, 
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initial  propellant  temperature  and  composition  of  the  propellant. 

(jo  2.4) 

Similar  work  has  been  conducted  at  Stanford  Research  Institute  '  ’  . 

It  was  observed  in  References  7  and  8  that  once  a  motor  is 
operating  in  an  unstable  manner,  its  behavior  is  independent  of  the 
process  which  initiated  the  instability.  Thus,  except  for  a  short 
initiation  transient,  artifical  triggering  produced  the  same  engine 

(7- 

response  found  in  the  spontaneous  cases.  In  the  above  experiments 
12) 

,  the  ultimate  test  for  stability  was  the  response  of  the  engine 

(25 ) 

flow  field  to  a  sharp  "pui.se"  v  '  created  by  an  explosive  charge 
placed  at  the  head  end  of  the  combustion  chamber.  Several  charges 
were  exploded  at  predetermined  times  during  the  motor  firing  while 
the  pressure-time  history  at  the  head  and  aft  ends  was  recorded. 

With  this  technique,  the  results  are  a  measure  of  the  nonlinear  sta¬ 
bility  of  the  motor;  i.e.,  does  a  large  amplitude  disturbance  decay 
and  the  flow  field  return  to  steady  state,  or  does  it  lead  to  an 
increasing  mean  chamber  pressure  and  burniiig  rate,  and  instability 

in  the  chamber?  In  a  recent  study  to  assess  the  effect  of  propellant 

(12) 

composition,  Roberts  and  Brownlee  attempted  to  correlate  the  onset 

of  unstable  operation  (eg.,  transition  pressure,  transition  burning 
rate)  with  propellant  ingredients  and  the  severity  of  the  disturbance. 
Quoting  from  Reference  (12): 

Clearly,  it  is  possible  to  investigate  a  particular  propellant 
system  and  to  establish  individual  trends  based  on  selected  changes 
to  the  formulation.  It  is  also  evident  that  such  trends  cannot  be 
expected  to  apply  to  propellants  containing  substantially  different 
binder  or  oxidizer  systems. 

The  concluding  statement  is :  "This  result  serves  to  underline  the  impor¬ 
tance  of  achieving  an  improved  understanding  of  the  physical  processes 
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governing  longitudinal  combustion  instability." 

The  purpose  of  this  report  is  to  develop  a  theoretical 
axial -mode  model  of  a  solid  propellant  rocket  engine  with  a  cylindri- 
cally  perforated  grain  which  can: 

(1)  predict  the  time -dependent  behavior  of  the  combustion 
chamber  and  nozzle  flow  field  after  the  introduction  of  an  arbitrary 
disturbance,  i.e.,  simulate  the  experimental  procedure  for  assessing 
nonlinear  stability,  and 

(2)  lead  to  a  better  understanding  of  the  complex  coupling 
between  the  chamber  flow  conditions  and  the  propellant  combustion 
process  in  longitudinal  combustion  instability. 

A  theoretical  description  of  the  nonlinear  behavior  of  axial 
instabilities  will  require  a  flow  field  model  with  certain  minimum 
capabilities.  The  need  for  a  solution  to  the  full  unsteady  conser¬ 
vation  equations  for  a  two-phase  flow  is  obvious.  Traveling  discon¬ 
tinuities  are  acceptable  solutions  to  these  nonlinear  equations  and 
hence  must  be  provided  for.  In  the  previously  mentioned  experiments, 
the  presence  of  traveling  shock  waves,  although  weak  (M  <  1.2),  is 
indicated  by  the  pressure- time  histories  and  confirmed  by  optical 
measurements  Since  the  growth  or  decay  of  disturbances  in  the 

combustion  chamber  is  often  decided  by  a  slight  imbalance  between  the 
overall  gain  and  loss  processes  of  the  system,  the  boundary  conditions 
imposed  upon  the  flow  field  model  become  very  important.  The  choked 
throat  of  the  nozzle  provides  the  physical  downstream  boundary  condition 
to  disturbances  in  the  chamber.  The  theoretical  model  must  account 

for  this  effect  without  imposing  any  restrictions  on  the  magnitude  of 
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the  mean  flow  Mach  number  or  the  amplitude  of  the  disturbance.  Equally 
important  is  the  solution  to  the  solid  propellant  combustion  problem 
which  has  boundary  conditions  dependent  on  the  chamber  flow  properties. 
This  solution  must  account  for  the  unsteadiness  of  the  burning  process 
without  making  assumptions  about  the  type  of  disturbance  in  the  chamber. 

Since  the  combustion  model  is  a  critical  part  of  an  analysis  of 
transient  behavior  inside  a  solid  propellant  rocket  motor,  it  warrants 
further  elaboration.  Regretably,  the  combustion  of  solid  propellants 
involves  complex  mechanisms  which  are  not  well  understood,  even  in 
steady  state.  Changes  in  composition  and  heterogeneity  have  substantial 
influence  on  the  measured  values  of  the  steady-state  surface  regression 
rate  as  a  function  of  pressure.  This  type  of  effect  cannot  be  predicted 
theoretically  with  current  models  of  the  combustion  process.  However, 
the  difficult  theoretical  problem  has  often  been  avoided  by  use  of  an 
empirical  expression  of  the  form, 

burning  rate  ot  pn 


or 


r  =  apn  (l.l) 

where  a  is  a  constant  and  n  is  called  the  pressure  index.  This  well- 
known  burning  rate  "law"  resulted  from  the  correlation  of  experimental 
steady  state  burning  rate  data  and  is  usually  applicable  within  a 
restricted  pressure  range.  Another  steady  state  expression, 
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(1.2) 


1 

r 


was  suggested  by  the  Granular  Diffusion  Flame  (GDF)  theory  for 

composite  propellants.  The  GDF  model  is  based  on  a  pre-raixed  kinetically 
controlled  reaction  at  low  pressure  (r  a  p)  and  a  diffusion  controlled 

2/0 

reaction  involving  "pockets  of  fuel  vapor"  at  high  pressure  (r  a  p  /J). 

f  17  ) 

With  suitable  values  of  a  and  b,  Equation  (1.2)  correlates  '  '  certain 
experimental  data  for  AP  -  hydrocarbon  binder  propellants  with  a  uni- 
modal  oxidizer  particle  distribution.  For  arbitrary  propellants  though, 
the  predictive  nature  of  both  Equation  (l.l)  and  (1.2)  must  be  viewed 
with  caution. 

If  the  burning  propellant  is  exposed  to  a  transient  flow,  as  in 
the  present  study,  the  use  of  any  steady  state  burning  rate  expression 
is  open  to  question.  This  difficulty  can  be  illustrated  with  a  simpli¬ 
fied  description  of  solid  propellant  combustion.  Consider  a  gaseous 

4 

flame  zone  adjacent  to  a  regressing  propellant  surface  (y  =  0)  in  a 
one-dimensional  situation  (Figure  l).  The  effect  of  heat  released  in 
the  flame  zone  and  at  the  surface  (if  surface  reactions  are  present) 
is  to  determine  the  temperature  gradient  in  the  solid  at  the  interface 
y  =  0.  With  this  boundary  condition,  the  time  dependent  energy  balance 
in  the  unburned  propellant  determines  the  temperature  at  the  interface. 
The  mass  flow  rate  from  the  surface  (or  surface  regression  rate)  is 
primarily  a  function  of  this  temperature.  Hence,  a  step  change  in 


Figure  1.  Typical  Solid  Propellant  Combustion  Model 

pressure  or  temperature  in  the  flame  zone  will  not  produce  an  instantane¬ 
ous  adjustment  in  surface  mass  flow  rate.  The  response  time  (time  delay) 
before  a  new  equilibrium  is  achieved  is  related  to  the  thermal  wave 
travel  time  in  the  unburned  propellant.  Use  of  a  steady  state  burning 
law  implies  a  constant  equilibrium  between  pressure  and  burning  rate, 
which  may  be  valid  only  in  restrictive  cases .  An  accurate  assessment 
of  whether  a  given  disturbance  in  the  combustion  chamber  will  grow  or 
decay  must  include  the  effect  of  the  combustion  response  time. 

Cheng  (-*-^"20)  hag  ana]_yZe(j  solid  propellant  combustion  with  a 

(21-23) 

time-lag  concept  inspired  by  Crocco's  '  work  in  liquid  propellant 

rocket  motors.  The  assumption  is  made  that  the  mass  flow  rate  from  the 
surface  is  proportional  to  the  pressure  at  a  previous  time,  i.e., 

m  (t)  ~  [p(t  -  T)]n  (1.3) 

o 
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The  time  lag  parameter,  t,  incorporates  the  effect  of  the  entire  com¬ 
bustion  process,  but  no  method  is  given  to  compute  its  value.  This 

feM  (25) 

problem  also  arises  in  the  works  of  Grad  *  and  Moore  and  Maslen  v  ' 

which  are  based  on  different  definitions  of  a  time  lag.  In  addition  to 
the  fact  that  "t"  will  be  strongly  dependent  on  the  propellant  pro¬ 
perties,  chamber  length,  type  of  disturbance  in  the  chamber,  flow  field 

composition  and  temperature,  it  may  also  change  sign.  More  recent 
(2&) 

analyses  '  '  of  the  solid  propellant  combustion  response  to  small 

amplitude  periodic  disturbances  in  pressure  have  shown  that  the  burning 
rate  can  lead  the  pressure  oscillation  at  low  frequencies.  Even  if  a 
straight-forward  method  to  compute  the  value  of  "t"  was  known,  the 
computational  difficulties  associated  with  storing  flow  fields  obtained 
at  previous  times  and  then  interpolating  in  time  and  space  to  determine 
the  proper  variable  at  a  retarded  time  would  be  immense.  Thus,  the 
time-lag  concept  will  not  be  considered  in  the  present  analysis. 

A  more  detailed  model  of  the  combustion  process  is  necessary  to 
circumvent  the  problem  of  an  unknown  time-lag  parameter.  The  gross 
features  of  the  process  can  be  represented  by  three  regions:  (l)  the 
cold  solid  region  in  which  a  thermal  wave  propagates  into  the  unburned 
propellant,  heating  the  material  to  the  surface  temperature,  (2)  an 
interface  region  between  the  gas  and  solid  phases  (often  treated  as 
infinitesimal  plane)  in  which  decomposition,  pyrolysis,  and  other  sur¬ 
face  reactions  may  occur,  and  (3)  a  gas-phase  reaction  zone  (flame). 

Inf  ) 

Many  theoretical  combustion  models  '  ' ,  based  on  this  interpertation 

of  the  combustion  process,  have  been  used  to  investigate  the  response 

of  the  burning  propellant  to  small  amplitude  periodic  pressure 
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disturbances.  Although  these  models  differ  as  to  the  treatment  of 
details,  they  share  the  following  simplifying  assumptions: 

(1)  a  one-dimensional  homogeneous  solid  phase 

(2)  simple  pyrolysis  of  solid  to  gas 

(3)  no  condensed  phase  reactions 

(4)  quasi-steady  gas  phase  flame  zone 

As  demonstrated  in  the  excellent  review  by  Culick.  these  common 

assumptions  are  responsible  for  each  analysis  arriving  at  the  same 

tu  *  /m 

two-parameter  form  of  the  combustion  response  function,  R  =  ,  as 

P'/p 

a  function  of  the  frequency  of  oscillation  of  the  disturbance.  This 
form  of  the  response  function  does  predict  that  the  surface  mass  flow 
rate  as  a  function  of  time  will  lead  a  periodic  pressure  disturbance 
at  low  frequencies  and  lag  at  high  frequencies.  It  also  predicts  a 
single  response  peak  or  maximum  at  a  non-zero  frequency  determined 
by  the  values  of  the  two  adjustable  parameters.  Qualitatively,  this 
behavior  corresponds  to  observed  experimental  trends,  but  quantitative 

(27) 

predictions  have  met  with  only  limited  success  .  Correlation  of 
all  experimental  data  will  undoubtedly  require  a  very  sophisticated 
theoretical  analysis  (eg.,  see  the  next  paragraph)  which  has  yet  to  be 
formulated.  However,  the  above  models  have  yielded  considerable  in¬ 
sight  into  the  transient  combustion  response  of  solid  propellants;  hence, 
their  basic  assumptions  will  be  adopted  in  the  present  study.  But  the 
response  function  derived  from  a  linearized  analysis  which  assumes  the 
existance  of  small  amplitude  periodic  disturbances  is  not  applicable 
when  arbitrary  disturbances  can  be  present  in  the  engine  flow  field. 

The  latter  situation  will  require  a  time-dependent  solution  of  the 

27 


non-linear  equations  which  follow  from  the  combustion  model. 

The  above  discussion  of  solid  propellant  combustion  has  ignored 
a  difficult  and  important  aspect  of  the  problem.  The  unsteady  propellant 
burning  rate  can  be  sensitive  to  velocity  fluctuations  parallel  to  the 
surface  as  well  as  pressure  fluctuations ''in  the  flame  zone.  The  former 


is  referred  to  as  velocity  coupling.  The  basis  for  velocity  coupling 
is  often  attributed  to  a  complexity  associated  with  the  combustion  of 
metal-loaded  solid  propellants  -  the  layer  of  molten  metal  on  the  pro¬ 
pellant  surface.  However,  this  is  not  to  imply  that  velocity  coupling 
must  be  absent  if  a  non-metalized  propellant  is  used.  Price  v  '  gives 
the  following  description  of  the  combustion  of  aluminized  propellants: 


Because  of  its  high  boiling  point,  aluminum  tends  to  accumulate  on 
the  burning  surface  in  a  condensed  phase,  protected  from  ignition 
by  a  thin  layer  of  aluminum  oxide.  The  aluminum  leaves  the  burning 
surface  primarily  as  agglomerates  in  the  30  -  300  y  diameter  range 
(details  depending  on  propellant  variables).  The  accumulation  pro¬ 
cess  on  the  surface  evidently  involves  sintering  of  individual 
particles  together,  probably  concurrently  with  attainment  of  the 
melting  point  of  the  metal.  This  is  followed  eventually  by  the 
complete  breakdown  of  the  oxide  "skin"  on  original  particles  and 
formation  of  reacting  agglomorate  droplets.  The  extent  of  reaction 
on  the  burning  surface  is  unknown  and  probably  dependent  on  ingredi¬ 
ent  and  formulation  variables.  Ignition  is  often  seen  to  occur 
concurrently  with  release  from  the  surface,  but  with  some  propellants 
the  luminous  metal  droplets  remain  for  some  time  on  the  surface,  and 
with  other  propellants  no  droplet  ignition  occurs  untill  after  separ¬ 
ation  from  the  surface. 


Regardless  of  where  ignition  of  aluminum  occurs,  most  of  its 
combustion  occurs  after  it  leaves  the  propellant  surface,  and 
proceeds  for  a  considerable  distance  from  the  surface. 

Certainly,  if  the  characteristic  tines  for  the  surface  accumulation  and 

the  metal  droplet  combustion  processes  are  the  same  order  of  magnitude 

as  the  travel  time  of  a  disturbance  in  the  combustion  chamber,  then 

these  processes  represent  additional  sources  of  energy  transfer  which 
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may  not  be  "in-phase"  with  the  flow  field  disturbance.  Hence,  they 

are  potential  driving  mechanisms  of  combustion  instability.  Since 

particle  damping  is  minimal  for  the  low  frequency  longitudinal  chamber 

modes,  it  is  likely  that  the  exothermic  metal  combustion  process  and 

the  interaction  between  velocity  fluctuations  and  the  molten  layer  on 

the  burning  propellant  surface  contribute  to  driving  axial  oscillations 

in  the  engine.  Recall  the  experimental  observation  that  increasing  the 

percent  of  metal  loading  in  the  propellant  composition  often  aggrevates 

the  longitudinal  combustion  instability  problem. 

Evidence  to  date  indicates  that  velocity  coupling  is  insensitive 

to  the  direction  of  the  velocity  fluctuations  in  the  chamber  flow  field 

Additionally,  its  contribution  to  the  propellant  response  is 

non-zero  only  when  the  magnitude  of  the  velocity  disturbance  exceeds 

(2) 

a  certain  threshold  value  .  Thus,  the  eoupling  is  distinctly  non¬ 
linear.  At  present,  no  satisfactory  theoretical  model  of  this  complex 

(2  29-31) 

interaction  is  available.  Besides  the  acoustic-type  analyses  ’ 

which  treat  the  effect  of  velocity  coupling,  two  computational  methods 

have  been  used  to  describe  this  phenomenon:  (l)  a  correlation  expres- 

sion  derived  originally  for  steady  erosive  burning,  and  (2)  use 

of  the  two-parameter  response  function  ^^’33)  £er;pyC(i  for  pressure 

(33  34) 

coupling,  but  modified  '  *  •  with  a  threshold  velocity  magnitude  term. 

Without  stronger  theoretical  and  physical  justifications,  the  applica¬ 
bility  of  these  computational  methods  are  open  to  question.  Therefore, 
the  description  of  the  unsteady  combustion  response  in  the  present 
study  will  be  limited  to  consideration  of  pressure  coupling  only, 

recognizing  that  the  results  obtained  for  axial-mode  instability  in 

29 


a  solid  propellant  rocket  motor  must  be  carefully  interperted. 

Previous  Work  in  Nonlinear  Instability 

The  observed  behavior  of  solid  rocket  motors  subjected  to  pulsing, 
and  the  fact  that  traveling  compression  wave  axial  instabilities  can 
exist  at  shock  strength  indicate  that  controlling  mechanisms  are 
definitely  nonlinear.  The  need  to  account  for  nonlinear  effects  has 
long  been  recognized,  but  the  solution  to  the  full  unsteady  conser¬ 
vation  equations  presents  a  formidable  prob?.em.  Furthermore,  certain 
details  of  the  combustion  process  are  not  well  enough  known  to  con¬ 
struct  a  realistic  mathematical  model.  As  a  result,  most  theoretical 
analyses  of  solid  rocket  instability  have  been  based  on  the  linearized 
equations  which  govern  the  behavior  of  small  amplitude  disturbances. 

Then,  the  influence  of  a  single  nonlinear  effect  is  studied  within  the 
framework  of  an  acoustic  chamber  analysis. 

Since  axial -mode  instabilities  involve  fluctuations  of  the  gas 

velocity  parallel  to  the  propellant  surface,  the  possible  contribution 

(29-31) 

of  velocity  coupling  received  early  attention  '  ^  .  In  Reference  29, 

the  combustion  response  due  to  velocity  coupling  is  assumed  proportional 
to  the  magnitude  of  the  difference  between  the  mean  flow  velocity  and 
the  instantaneous  velocity.  This  rectification  effect  introduces  a 
nonlinearity  proportional  to  the  first  order  velocity  perturbation, 
whereas  the  nonlinearity  associate;,  with  the  pressure  is  proportional 
to  a  squared  perturbation.  Use  of  a  linearized  analysis  in  the  chamber 
is  justified  on  this  basis.  The  possibility  of  attenuation  or  ampli¬ 
fication  of  disturbances  adjacent  to  the  burning  propellant  is  discussed 
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in  terras  of  the  values  of  the  pressure-coupled  response  function  and 

a  "frequency-dependent  erosion  constant",  both  of  which  are  unknown. 

(2) 

Price  and  Dchority  have  examined  a  similar  problem  but  introduced 
the  further  nonlinear  effect  of  a  threshold  velocity,  i.e.,  the  com¬ 
bustion  response  is  non-zero  only  when  the  stream  velocity  exceeds  the 
threshold  value  (analogous  to  erosive  burning).  This  analysis  gives 
considerable  insight  into  how  velocity  coupling  should  affect  oscilla¬ 
tions  in  the  chamber,  but  it  is  not  detailed  enough  to  predict  the 

(qq) 

growth  of  a  disturbance  to  a  limiting  amplitude,  etc..  Culick 
looked  at  the  stability  of  longitudinal  oscillations  with  a  combined 
wave  equation  valid  in  the  chamber.  The  nonlinear  velocity  coupling 
phenomenon,  described  in  a  similar  manner  as  in  Reference  (2),  enters 
as  a  boundary  term  to  the  linear  analysis. 

An  entirely  different  type  of  analysis  was  undertaken  by  Marxman 

('ll, ) 

v  to  explain  the  existence  of  traveling  shock  wave  instabilities. 
Without  combustion  driving,  a  traveling  shock  wave  would  quickly  be 
overtaken  by  expansion  waves  from  the  high  pressure  side.  Thus,  it 
is  theorized  that  as  the  shock  wave  passes  over  the  propellant,  the 
increase  in  surface  mass  flow  rate  must  be  sufficient  to  keep  the  down¬ 
stream  flow  sonic  relative  to  the  wave  (similar  to  a  detonation  wave). 
Using  a  resjjonse  function  for  pressure  coupling,  the  minimum  value  of 
combustion  response  required  to  sustain  the  traveling  shock  is  estimated 
from  a  local  analysis.  This  estimate  is  in  general  agreement  with  the 
results  found  in  the  present  study.  However,  Marxman1  s  work  does  not 
account  for  mean  flow  effects,  the  nozzle  loss,  or  transient  behavior. 
The  first  analysis  to  include  coupled  nonlinear  effects  without 
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recourse  to  order  of  magnitude  estimates  for  eliminating  various  terms 

(34) 

in  the  equations  is  due  to  Levine  and  Culick  '  ' .  This  investigation 

is  based  on  a  numerical  solution  of  the  full  time-dependent  conser¬ 
vation  equations  for  a  two-phase  flow,  and  hence  represents  a  new 
approach  in  the  solution  of  instability  problems.  The  quasi-one- 
dimensional  model  in  Reference  (34)  invokes  the  standard  assumptions 
of  an  inviscid,  perfect  gas  and  inert,  mono-disperse,  constant  diameter, 
spherical  particles .  The  solution  for  the  chamber  flow  is  formulated 
as  an  initial-value  problem  and  solved  with  the  method-of-characteristics . 
The  chamber  boundary  condition  at  the  nozzle  entrance  plane  is  derived 
from  an  isentropic  relationship  between  the  flow  properties  at  the  nozzle 
entrance  and  those  at  the  sonic  throat.  The  combustion  response  due  to 

(26) 

pressure  coupling  is  based  on  the  well-known  two-parameter  '  ’  response 

function  obtained  in  linearized  combustion  analyses  which  assume  small 

amplitude  periodic  motion  in  the  gas  phase.  The  desired  combustion 

response  is  obtained  in  real  time  by  evaluating  the  inverse  Laplace 

transform  of  the  given  function  of  frequency.  Velocity  coupling  is 

(33) 

treated  in  an  analogous  manner  v  .  The  same  response  function  but 
with  different  adjustable  parameters  and  a  multiplicative  coefficient 
to  account  for  the  threshold -magnitude -effect  is  used  to  predict  the 
response  to  velocity  fluctuations.  The  total  combustion  response  is 
taken  as  the  sum  of  the  two.  However,  since  the  two-parameter  response 
function  was  derived  for  small  amplitude  periodic  disturbances,  its 
use  in  a  combustion  chamber  with  arbitrary  amplitude  transient  dis¬ 
turbances  raises  the  following  questions : 

(l)  How  can  this  model  determine  the  transient  response  of  the 
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propellant  before  a  definite  frequency  of  oscillation  has  been  estab¬ 
lished  in  the  chamber,  and 

(2)  How  can  a  linear  model  determine  the  combustion  response 
to  finite  amplitude  disturbances? 

Statement  of  the  Problem 

It  is  the  objective  of  this  investigation  to  develop  a  theoretical 
one-dimensional  model  for  solid  propellant  rocket  engines  with  a  cylin- 
drically  perforated  grain  which  can  predict  the  time-dependent  behavior 
of  the  combustion  chamber  and  nozzle  flow  field  after  the  introduction 
of  an  arbitrary  disturbance.  The  solution  will  include  the  unsteady 
flow  in  the  choked  nozzle,  eliminating  the  need  for  an  artificial  nozzle 
boundary  condition.  The  combustion  of  the  solid  propellant  will  be 
described  by  a  time-dependent  solution  obtained  simultaneously  with 
the  chamber  flow  field.  No  restrictions  vail  be  placed  on  the  amplitude 
or  type  of  disturbance  in  the  engine.  Shock  waves  resulting  from  com¬ 
bustion  driving  or  as  an  initial  disturbance  simulating  the  experi¬ 
mental  pulsing  technique,  will  be  acceptable  solutions  to  the  theoretical 
model . 

In  what  follows,  Chapter  II  describes  the  development  of  the 
solutions  to  the  nonlinear-  time-dependent  conservation  equations  for 
a  two-phase  flow  in  the  combustion  chamber  and  nozzle .  Chapter  III 
describes  the  analysis  and  solution  of  the  time-dependent  solid  pro¬ 
pellant  combustion  process  which  will  be  used  to  predict  the  combustion 
response  to  arbitrary  disturbances  in  the  chamber.  Chapter  IV  discusses 
the  results  obtained  for  the  transient  response  of  the  engine  after 


introduction  of  various  types  of  disturbances.  Chapter  V  summarizes 
the  results  obtained  in  this  investigation. 


! 

I 


34 


CHAPTER  II 


ANALYSIS  OF  THE  COMBUSTION  CHAMBER  AND 
NOZZLE  FLOW  FIELD 

Discussion  of  Assumptions 

Solid  propellant  rocket  motors  exhibiting  axial  mode  instability- 
are  usually  characterized  by  a  ratio  of  chamber  length  to  chamber  radius 
much  greater  than  unity.  Considering  the  high  turbulence  level  present, 
radial  property  variations  are  sufficiently  small  that  a  one-dimensional 
formulation  provides  an  adequate  description  of  the  flow  field.  To 
avoid  placing  restrictions  on  the  type  of  transient,  behavior  which  may 
be  encountered  in  the  rocket  motor,  the  present  investigation  will 
cons ider, 

(a)  the  full  non-linear,  time -dependent  conservation  equations 

(b)  the  flow  in  the  choked  nozzle 

(c)  the  transient  behavior  of  the  combustion  process  coupled 
to  the  chamber  flow  field  (see  Chapter  III). 

The  first  requirement  means  that  the  analysis  must  be  capable  of  handling 
traveling  shock  waves,  if  they  appear  in  the  flow  field. 

With  these  requirements,  consider  an  axisymmetric,  one -dimensional 
combustion  chamber  and  nozzle  system  with  a  cylindrically  perforated 
solid  propellant  grain  (see  Figure  2).  The  solid  propellant  is  assumed 
to  be  the  metal  loaded  type  (typically,  aluminum)  which  deposits  solid 
particles  (eg.,  Al^O^)  in  the  flow  when  combustion  is  completed.  Hence, 
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mass  flow  rate  of  gas 
m9“  unit  perimeter  area 

mass  flow  rate  of  particles 

mp=  — - : — - - c - 

unit  perimeter  area 

h*  =  enthalpy  of  combustion 


Figure  2.  Solid  Propellant  Rocket  Engine  Model 


the  flow  in  the  combustion  chamber  and  nozzle  is  two  phase.  The  present 
analysis  is  based  on  several  standard  assumptions: 

(1)  The  gas  medium  is  thermally  and  calorically  perfect.  Implied 
in  this  statement  is  the  assumption  that  all  combustion  is 
complete  when  the  gases  have  left  the  flame  zone  and  entered 
the  combustor  control  volume. 

(2)  The  gas  medium  is  inviscid  and  non-heat-conducting.  However, 
the  gas -particle  interaction  can  involve  the  exchange  of 
momentum  and  energy.  The  heat  transfer  process  at  the 
burning  propellant  surface  located  outside  the  combustor 
control  volume  perimeter  boundary  is  dealt  with  separately 
in  Chapter  III. 
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(3)  The  solid  particles  suspended  ii  the  gas  have  negligible 
volume  and  exert  no  pressure. 

(4)  The  actual  distribution  of  particle  sizes  can  be  replaced 

with  a  system  of  non-reacting  spherical  particles  having 
one  average  diameter  constant  in  time.  As  Price  has 

emphasized,  the  solid  particles  are  not  produced  directly 
in  the  flame  zone.  Since  the  ignition  temperature  is  much 
higher  than  the  melting  point  temperature  for  metals,  a 
liquid  layer  forms  on  the  propellant  surface.  Complicated 
condensed  phase  reactions  occur  between  the  gas  and  large 
agglomerates  of  molten  metal  from  the  surface.  These 
reactions  may  take  place  well  outside  of  the  flame  as  the 
rnetal  droplets  are  entrained  in  the  gas  flow.  When  com¬ 
bustion  is  complete,  solid  oxide  particles  remain  suspended 
in  the  gas.  In  general,  the  resultant  distribution  of 
particle  diameters  is  bimodol;  that  is,  particles  in  the 

0  -  (smoke)  range  and  the  larger  5  -  20p,  range.  The 
approximate  nature  of  assumption  (4)  is  fully  recognized, 
but  further  modeling  of  this  complex  process  is  beyond 
the  scope  of  the  present  analysis. 

(5)  The  drag  force  between  solid  particles  and  gas  can  be 

represented  by  a  simple  drag  law  as  a  function  of  Reynolds 

number.  In  view  of  the  discussion  of  assumption  (4),  use 

of  a  drag  law  for  spherical,  non-reacting  particles  in  the 

present  problem  is  also  open  to  question,  however,  the 

well-known  Stokes  Flow  Drag  Law  is  employed  in  this  study, 
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recognizing  that  the  upper  Reynolds  number  limit  may  be 

violated  in  certain  instances.  The  available  higher  order 

(34  35) 

corrections  ’  '  could  easily  be  introduced  in  future 

work. 

(6)  The  time  rate  of  change  of  thermal  energy  stored  by  the 
particles  is  negligible  compared  to  the  time  rate  of  change 
of  kinetic  energy.  Thus,  particles  are  assumed  to  leave 
the  nozzle  with  whatever  thermal  energy  they  acquired  in 
the  flame  zone,  which  is  valid  for  a  vanishingly  small 
heat  transfer  coefficient.  Reference  (34)  has  not  made 
this  assumption  and  has  included  the  energy  transfer  between 
gas  and  particles.  It  appears,  however,  that  this  effect 

is  small  and  its  treatment  could  be  delayed  until  a  more 
sophisticated  combustion  analysis  is  attempted. 

(7)  The  mass  flow  rates  of  gas  and  particles  from  the  burning 
propellant  surface  enter  the  perimeter  boundary  of  the 
combustor  control  volume  with  negligible  velocity.  Thus, 
the  mass  addition  at  the  perimeter  boundary  contributes 
zero  axial  momentum  and  zero  kinetic  energy  to  the  control 
volume.  Part  of  the  energy  in  the  control  volume  must  be 
used  to  accelerate  the  entering  mass  flows  of  gas  and 
particles  to  their  respective  local  velocities .  This 
results  in  an  additional  drag  term  in  the  momentum  equations 
and  a  dissipation  term  in  the  energy  equation,  which  approxi¬ 
mates  the  effect  of  a  wall  boundary  layer. 

(8)  The  increase  of  chamber  volume  due  to  regression  of  the 
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propellant  surface  is  assumed  negligible.  For  a  linear 
regression  rate  of  one  cm/sec  and  an  initial  chamber  radius 
of  one  foot,  the  volume  change  after  one  second  is  about 
six  percent.  Since  wave  travel  times  are  on  the  order  of 
milliseconds,  this  effect  can  be  ignored.  The  exception 
might  be  ignition  transient  phenomena  if  several  seconds 
are  involved. 

The  one -dimensional  assumption  needs  further  comment.  The  con¬ 
trol  volume  shown  in  Figure  2  illustrates  the  inherent  difficulty  with 
a  one -dimensional  analysis  as  applied  to  this  problem,  'the  mass  flow 
of  gas  and  particles  leaving  the  propellant  surface  enters  perpendicular 
to  the  axial  direction.  In  the  present  analysis,  the  integral  of  the 
ma  s  flow  rate  over  the  perimeter  boundary  of  the  control  volume  becomes 
a  mass  source  within  the  control  volume.  However,  the  exact  process 
by  which  the  mass  flows  of  gas  and  particles  are  accelerated  to  their 
respective  axial  velocities  within  the  control  volume  is  obscured  by 

the  integral  approach.  Although  the  process  is  left  unspecified,  its 

(34) 

effect  is  accounted  for  in  the  analysis.  Levine  and  Culick  '  '  employ 

the  same  approximation.  This  avoids  a  difficult  two-dimensional  analy¬ 
tical  problem  but  the  higher  order  formulation  will  eventually  be 
required  to  model  the  details  of  the  interaction  of  the  "core"  flow 
with  the  combustion  process  at  the  boundary. 

Method  of  Solution 

With  the  assumptions  of  the  previous  section,  the  time-dependent 
convervation  equations  applied  to  the  present  problem  lead  to  a  hyper- 
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bolic  initial -value  problem.  Solutions  to  this  system  of  equations  can 
include  traveling  discontinuities.  When  developing  a  method  of  solution, 
an  important  consideration  is  the  accurate  tracking  of  these  discon¬ 
tinuities. 

In  computational  fluid  mechanics,  two  basic  methods  of  employing 
the  conservation  equations  have  been  advanced.  One  widely-used  technique 
(36-38)^  referreii  to  as  the  finite-element  method  or  conservative  form, 
is  one  in  which  mass,  momentum  and  energy  are  conserved  strictly  for 
a  finite-size  control  volume.  An  advantage  '  of  this  method  is 
attributed  to  a  cancellation  of  numerical  error  on  interior  boundaries 
when  elements  are  summed  over  a  finite  domain.  Since  this  formulation 
reduces  to  the  proper  jump  conditions  as  the  computational  element 
shrinks  to  zero  thickness,  it  has  been  claimed  u  '  that  discontinuous 

behavior  can  be  computed  without  special  computational  procedures. 

(39) 

Using  MacCormack's  '  predictor-corrector  integration  scheme,  Kutler 

(oQ) 

and  Lomax  '  have  computed  steady  flow  fields  (which  include  shock 
waves)  over  various  blunt  body  and  wing  configurations.  Their  method 
is  said  to  be  "shock  capturing",  i.e.,  shock  waves  are  smeared  over 
three  or  more  mesh  points  but  automatically  forced  to  the  proper 
location.  A  second  method  is  to  solve  the  conservation  equations  in 
partial  differential  form  with  a  finite  difference  technique.  Then 
the  information  contained  in  the  associated  characteristic  equations 
can  be  used  to  track  discontinuity  5  precisely.  Moretti  has 

used  this  method  for  shock  tube  problems  with  multiple  shock  waves  and 
contact  discontinuities . 

In  Reference  (4l),  the  familiar  problem  of  a  piston  accelerating 
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into  a  stagnate  gas  is  used  to  compare  the  predictions  of  many  well- 
known  numerical  integration  schemes,  employing  both  the  conservative 
and  finite  difference  form  of  the  conservation  equations.  For  the 
time  period  before  shock  wave  formation,  second-order  integration 
schemes  with  the  conservative  form  of  the  equations  demonstrate  no  more 
accuracy  than  with  the  finite  difference  form.  All  methods  yield 
acceptable  results.  However,  if  the  computation  is  continued  beyond 
the  formation  of  the  shock  wave,  numerical  waggles  appear  in  the  flow 
field  with  all  schemes.  A  numerical  wiggle  is  defined  here  as  a  "saw¬ 
tooth"  oscillation  (high,  low,  high,  low,  .  .  . )  of  the  magnitude  of 
a  flow  field  variable  from  one  mesh  point  to  the  next.  Doubling  the 
number  of  mesh  points  usually  doubles  the  number  of  wiggles  since  no 
physical  process  is  represented.  Although  MacCormack's  scheme  with 

the  conservative  form  of  the  equations  (i.e.,  the  shock  capturing  method) 
exhibits  the  minimum  amount  of  this  behavior,  pronounced  wiggles  are 
evident  on  the  low  pressure  side  of  the  shock  wave  even  with  360 

mesh  points  between  the  piston  and  the  undisturbed  portion  of  the  gas. 

By  contrast,  Moretti  solved  the  same  problem  with  20  mesh  points 

using  the  finite  difference  form  of  a  second-order  Taylor  series 
expansion  of  the  conservation  equations  and  supplementary  information 
from  the  method-of- characteristics,  and  predicted  the  behavior  of  the 
discrete  shock  wave  with  no  wiggles.  It  is  concluded  that  integrating 
continuous  differential  equations  across  a  discontinuity  will  produce 
wiggles . 

The  properties  demonstrated  by  Moretti's  computational  scheme 
are  required  in  the  present  transient  analysis  where  the  location  of 


a  shock  wave  in  the  rocket  engine  flow  field  is  important  in  deter¬ 
mining  the  instantaneous  combustion  response.  Consequently,  Moretti's 
(In  42) 

v  5  '  technique  is  adopted  in  this  study.  The  method  is  based  on 

three  important  postulates  which  are  paraphrased  here: 

(1)  If  a  region  of  continuous  flow  is  defined  by  continuity 
of  second  derivatives,  then  all  integration  schemes  of  second-order 
accuracy  will  predict  acceptable  results  in  a  continuous  region.  Under 
these  conditions,  the  conservative  form  of  the  conservation  equations 
is  no  more  accurate  than  the  finite  difference  form. 

(2)  Treat  discontinuities  as  discontinuities;  i.e.,  do  not  4 
integrate  continuous  differential  equations  across  steep-fronted  waves. 

(3)  Do  not  violate  the  conservation  equations  when  evaluating 
flow  field  variables  on  a  boundary. 

These  principles  will  be  implemented  in  the  present  numerical  inte¬ 
gration  scheme  which  will  be  outlined  in  the  next  two  sections .  In 
a  brief  summary,  the  conservation  equations  will  be  written  for  a 
general  stream-tube  control  volume  and  then  reduced  to  partial  differ¬ 
ential  form  by  shrinking  the  control  volume  thickness  to  zero.  These 
equations  will  then  be  transformed  to  a  floating  computational  coordinate 
system  which  treats  discontinuities  as  moving  internal  boundaries.  For 
this  formulation,  a  hybrid  integration  scheme  arises  naturally.  All 
regions  of  continuous  flow  between  two  boundaries  will  be  integrated 
forward  in  time  with  a  second-order  Taylor  series  expansion  of  the 
differential  equations.  Values  of  flow  field  variables  on  all  boundaries 
as  well  as  the  movement  of  the  internal  boundaries  will  be  obtained 

with  the  supplemental  information  from  the  method  of  characteristics. 
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All  necessary  equations  are  derived  in  the  two  sections  which  follow. 

By  way  of  comparison,  the  instability  analysis  of  Levine  and 

(34) 

Culick  uses  the  method-of-characteristics  to  obtain  the  entire 

solution.  Thus,  the  numerical  accuracy  of  the  present  analysis 
should  be  equal  to  that  of  Reference  (3*0*  However,  in  a  flow  field 
with  traveling  regions  of  compression  and  expansion,  the  natural 
characteristic  mesh  system  becomes  very  irregular  and  can  require 
a  lot  of  time-consuming  interpolation.  This  is  avoided  with  Moretti's 
method. 


Derivation  of  Equations 

Based  on  the  assumptions  discussed  in  the  first  section,  the 
conservation  of  mass,  momentum  and  energy  for  the  general  control 
volume  shown  in  Figure  2  can  be  written  as  follows : 

Conservation  of  Mass  (gas  phase) 


1 1 1 J p4V  + 1 J  ^  ^  ■  j  J 


Th  dA 
g  P 


(2.1) 


c.v. 


C  •  o  • 


c.s.p. 


Conservation  of  Mass  (solid  particles) 


jb_  P  P 
fit  J  J  J 


r  r 

p  dV  +  j 

p  J 


P  u 
KP  P 


dA 


■  j: 


til  dA 
P  P 


(2.2) 


c.v. 


c.s. 


c.s.p. 
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Conservation  of  Momentum  (gas  phase  and  solid  particles) 


£  J 1 1 (pu  +  vV dV  + 1 J 


(puu 


p  u  u  ) 
p  p  p 


dA 


=  Tf  (2.3) 


c.v. 


c.s. 


Conservation  of  Momentum  (solid  particles  alone) 


J  J  J  ppup  av  +  J"  J*  Wp 

c.v.  c.s. 


J  J  J  KPpfr'V  dV  (2-4) 


Conservation  of  Energy  (gas  phase  and  solid  particles) 


» J J 1  {p  [  +  t  ] +  pP  -£■ } dV 

c.v. 


+  JJ  {p[^r  +  T]“+pp"f'“p}’dS 

c.s. 


-!! 


-7®  dA 
Y-l  g  P 


(2.5) 


C  •  S  iP  • 
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Equation  of  State  (gas  phase) 


p  =  pe  =  ph 


(2.6) 


where  all  variables  have  been  non-cbimensionalized  with  respect  to  the 
following  reference  state, 


* 

Pr 


p*  R 
Kr  gas 


T* 


Y  R  T* 
T  gas  r 


T*  =  h*/c  =  e*/c  (2.7) 

r  r  p  r'  v 

and  the  area  A*  ,  and  length  L*.  The  right-hand  side  of  Equation  (2.3) 
ch 

represents  the  reaction  of  the  surroundings  on  the  control  surface  due 
to  internal  pressure  forces.  The  right-hand  side  of  Equation  (2.4) 
represents  the  volumetric  drag  due  to  the  presence  of  solid  particles 
traveling  with  velocity  u^.  The  parameter  K, 


K 


(2.8) 
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results  from  the  application  of  the  Stokes  Flow  Drag  Law. 


By  shrinking  the  control  volume  to  zero  thickness,  Equations  (2.1) 
through  (2.5)  can  be  written  in  the  partial  differential  form  needed 
for  a  Taylor  series  expansion.  After  some  manipulation,  these  equations 
become, 

Conservation  of  Mass  (gas  phase) 


2m 


+  0  &  +  u  &  =  r,  ■  -*  - 


to 


{(X 


to 


pu 


R 


dlnA 

dx 


(2.9) 


Conservation  of  Mass  (solid  particles) 


to  to  6p_  2m 

bt  Kp  to  p  to  2  ^ 


-  P^u. 


dlnA 


p  p  dx 


(2.10) 


Conservation  of  Momentum  (gas  phase) 


to  to  yp 


^  H  -  b3  s  ■  + 1  *gu  p  (2-u) 


Conservation  of  Momentum  (solid  particles) 
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+  |  (ihgu2  +  mpup2)  ]  +  K  pp(u-up)2  }  (2.13) 

where 

h  =  p/p 

and 

S  =  lnp  -  ylnp  (2.l4) 

In  Equation  (2.l4),  the  reference  state  is  assigned  a  zero  entropy  value. 

These  equations  are  employed  in  a  computational  coordinate  system 
which  contains  two  fixed  boundaries  (the  head  end  of  the  combustion 
chamber,  and  the  nozzle  supersonic  exit  plane)  and  any  number  (which 
may  be  zero)  of  internal  moving  boundaries  as  shown  in  Figure  3.  A 


shock  wave  discontinuity  is  treated  as  a  moving  boundary,  and  the  flow 
field  in  a  region  between  any  two  boundaries  is  assumed  continuous.  In 
the  general  case,  the  two  boundaries  of  a  flow  field  region  may  move 
independently  of  each  other.  To  avoid  having  to  compute  the  movement 
of  a  boundary  over  a  system  of  fixed  mesh  points  and  to  simplify  the 
calculation  procedure,  the  independent  variables  are  normalized  for 
one  region.  These  computational  coordinates  will  then  "float"  with 
the  movement  of  the  boundaries.  Consider  the  following  transformation: 

(1)  First,  let  §  =  |(x)  be  a  stretching  function  which  allows 
equally  spaced  grid  points  in  f  to  represent  a  non-uniform  distribution 
in  x.  With  the  proper  weighting,  this  function  permits  a  minimum  total 
number  of  mesh  points  to  accurately  represent  regions  where  flow  field 
properties  are  changing  rapidly  (e.g.,  near  the  choked  throat). 

(2)  Then,  in  §-space,  let  the  instantaneous  location  of  the 
right  hand  boundary  be  given  by  c(t),  and  the  left  hand  boundary  by 
b(t).  For  a  typical  region  of  continuous  flow,  define 


(2.15) 


T  s  t  - 


where  X  =  0  and  X  =  1  locate  the  left  and  right  hand  boundaries  res¬ 
pectively  (see  Figure  3). 
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Figure  3.  Computational  Coordinate  System 


The  partial  derivatives  transform  as 


JL  „  ± 

(2.16) 


b_  ±  J_  ± 

bt  “  e>T  c-b  &X 


r,  = 

11  c-b 


where 


e  -  (x-D  ^  b(t)  -  x  £  c(t) 


By  letting, 


Us,nu  +  T^bT 


(2.17) 


U 

P 


+  (c-b) 


(2.18) 


then  Equations  (2.9)  through  (2.13)  transform  to, 


ilfi+U-^+np—  =  R 

bT  &X  l|P  bX  1 


bp„  bp„  bu 

- 2  +  U  - 2  +  ftp  __£  =  R 

bT  p  bX  l|Pp  bX  2 


(2.19) 


(2.20) 


1  it  bu  *Ti  bp  _  tj 

w  +  u  bx  +  YP  6X  -  r3 


(2.21) 


bu  bu 

—2  +  u  — 2  =  R. 

bT  p  bX  5 


(2.22) 


SO 


i 


6S  +  _  R 

6T  U  6X  R4 


(2.23) 


In  the  special  case  when  all  internal  boundaries  are  absent,  c  -  b  =  1, 
and  the  stretching  function  is  linear,  7]  =  d|/dx  =  1,  then  U  =  u,  U  = 
Up  (since  dc/dt  =  db/dt  =  0  for  the  stationary  boundaries)  and  the 
above  set  of  equations  reduces  to  the  original  form,  Equation  (2.9) 
through  Equation  (2.13). 

Forward  time  integration  of  Equations  (2.19)  through  (2.23)  in 
a  continuous  region  is  accomplished  with  a  finite-difference  approxi¬ 
mation  to  a  Taylor  series  expansion,  truncated  after  the  second-order 
terms.  This  insures  an  accurate  representation  of  the  non-homogeneous 

parts  of  the  equations.  It  is  not  clear  that  a  straight-forward  appli- 

(39) 

cation  of  well-known  schemes  such  as  MacCormack's  or  the  two-step 

(43) 

Lax-Wendroff  '  would  provide  the  proper  accuracy  for  these  non- 
homogeneous  terms.  The  Taylor  series  expansion  for  the  dependent 
variables  at  time  T  =  Tq  +  AT  (denoted  by  subscript  l)  in  terms  of 
the  known  information  at  Tq  (denoted  by  subscript  o)  is 


P1  ”  po  +  pT^o  AT  +  PTT^ 


tiP 


0  2 


(2.24) 


AT 

p  =  p  +  p  1  AT  +  p  I 

»i  po  V0  V°  2 


(2.25) 
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1 


U, 


=  +  UTI0  AT  +  u^l 


TT  o 


Al^ 

2 


(2.26) 


+  u 

t» 


PT'0 


AT  +  u 

T» 


PTT'o 


(2.27) 


+ 


^T^o  ^  +  ^rT'rr 


TT 1  o 


(2.28) 


The  first-order  time  derivatives  follow  directly  from  Equations  (2.19) 
through  (2.23)  with  a  centered  differencing  method  for  spatial  derivatives. 
Thus  for  the  mesh  system, 


m-1  m  m+1 


Equation  (2.19)  at  point  "m"  becomes 


The  second  order  time  derivatives  are  obtained  from  Equation  (2.19) 
through  (2.23)  by  an  additional  time  differentiation.  From  Equation 
(2.19)  it  follows  that 
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-  IVX + 


u%  + 


PUTX? 


+ 


(2.30) 


where  a  terra  like  can  be  evaluated  from  Equation  (2.21)  as 


(u^x  =  ^  +  UuXX  +  vp  [pxx  -  PXpX/p]  +  \px/Yp}  (2.31) 

Second  order  spatial  derivatives,  such  as  Uyy,  are  simply 


u 


u, 


m+1 


+  u 


m-1 


-  2u 


m 


XX 


AX 


(2.32) 


By  continuing  this  development,  the  right-hand  sides  of  Equations  (2.24) 

through  (2.28)  can  be  expressed  entirely  in  terms  of  the  known  solution 

at  time  Tq.  Given  the  proper  AT,  these  equations  determine  the  new 

flow  field  variables  at  time  T  +  AT  for  all  interior  mesh  points, 

o 

as  shown  in  Figure  4. 

However,  the  Taylor  scries  expansion  does  not  yield  values  of 
flow  field  variables  on  a  boundary.  Determining  these  boundary  values 
is  the  subject  of  the  next  section. 
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Figure  4.  Region  of  Applicability  of  Taylor  Series  Integration 


Boundary  Conditions  and  Shock  Waves 
The  information  contained  in  the  system  of  Equations  (2.19) 
through  (2.23)  when  written  in  characteristic  form  will  be  used  to 
track  the  moving  boundaries  and  establish  values  of  flow  field  variables 
on  fixed  boundaries.  If  the  pressure  gradient  term  in  Equation  (2. 21) 
is  replaced  by 


&  = 
6X 


P 


(2.33) 


then  the  system  becomes  five  equations  in  the  five  dependent  variables 
p,  Pp>  u,  ^  and  S.  The  five  characteristic  directions  in  the  trans¬ 
formed  coordinate  system  are  found  to  be, 
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OK 

dT 


U 

U  -i-  7|a 
U  -  T)a 

Up  (double  root ) 


(2.34) 


where  a  =  p/p.  If  the  equation  describing  energy  transfer  between 


gas  and  particles  had  been  included  in  the  system,  then  would  become 
a  triple  root.  The  general  compatibility  relation  is  given  by. 


/  6X 
\  dT 


du 

dT 


(2.35) 


For  the  individual  roots,  Equation  (2.35)  yields: 
clX 

(a)  along  =  U  (gas  path  stream  line) 


dS  =  dT 


(2.36) 


which  is  the  energy  equation. 

clX 

(b)  along  —  =  U  ±  T]a  (left  and  right -running  characteristics) 
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R  R 

(^)dp±du={±R3  +  a(-i+^)}dT  (2.37) 


For  the  special  case  of  no  heat  or  mass  addition  from  combustion,  single 
phase  flow,  and  stationary  coordinates,  Equations  (2.36)  and  (2.37) 
reduce  to  the  familar  gas  dynamic  relations, 


dS 


0 


along 


dx 

dt 


u 


(  f  )  ±  Y  (  IT  )  =  0  along  dt  =  u  *  a  (2‘38) 

A  peculiar  difficulty  arises  along  the  characteristic  direction 
dX/dT  =  U  (dov  root).  The  general  compatibility  relation  (i.e.. 
Equation  2.35)  yields  only  the  trivial  result,  0=0.  Normally,  this 
indicates  that  two  equations  in  the  system  are  already  in  characteristic 
form.  This  is  true  for  the  particle  momentum  equation  (i.e.,  Equation 
2.22)  which  can  be  written  as, 


du  du  du 

+  u  —I £  -  — E  = 

dT  p  5X  dT 


Kr 


(2.39) 
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along  cLX/dT  =  U^.  However,  the  particle  continuity  equation  (i.e.. 
Equation  2.20')  cannot  be  written  in  characteristic  fora.  The  contra- 
diciton  can  be  traced  to  the  lack  of  an  equation  of  state  for  the 
solid  particle  flow.  With  no  mechanism  relating  particle  density  to 
pressure,  and  hence  no  particle  speed  of  sound,  the  particle  momentum 
equation  is  not  coupled  to  the  particle  continuity  equation.  The 
present  analysis  will  interpret  this  to  mean  that  Equation  (2. 20)  in 
the  form. 


dp  du 

dT^=  r2  ■  %  ax2,  (2-4o) 

can  be  integrated  along  <]X/dT  r=  U  after  an  integration  of  Equation 
(2.39).  The  results  from  Equation  (2.39)  can  be  used  to  estimate  an 
average  value  of  Su^/SX.  Reference  (3*0  handled  this  problem  in  a 
similar  manner. 

In  the  X-T  computational  space  of  the  previous  section,  the 
flow  field  representing  the  combustion  chamber  and  nozzle  has  two 
stationary  boundaries  and  any  number  of  internal  moving  boundaries . 

As  an  illustration,  consider  a  flow  field  with  one  right  running  shock 
wave  (Figure  5).  The  head  end  of  the  combustion  chamber  is  X  =  0  in 
region  [l],  the  right-running  shock  wave  is  the  boundary  between  .1] 
and  [2],  and  the  exit  plane  of  the  nozzle  is  X  =  1  in  region  [2], 

Two  boundary  conditions  are  known  at  the  head  end  of  the  com¬ 


bustion  chamber,  namely 
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Figure  5.  Computational  Coordinates  for  Example  Problem 
With  One  Right -Running  Shock  Wave 


U>X=0  =  Up'X=0  =  0 


(2.41) 


However,  the  numerical  computation  also  requires  values  of  p,  p^,  and 
S  at  this  location  at  every  time  step.  These  are  determined  with  the 
method-of -characteristics.  From  the  T  -  X  diagram  shown  in  Figure  6, 

T  t© 

_  ,  left  running  characteristic 


T  +  AT 
o 


head  end  of  the 
combustion  chamber 


/  i£=  u 


Region  [l] 


Figure  6.  Method  Used  to  Compute  Boundary  Values 
at  Head  End  of  Combustion  Chamber 


.i 


it  can  be  seen  that  the  left-running  characteristic  reaches  the  new 
wall  point  (2)  from  some  location  along  the  known  datum  plane  T  =  T  . 

The  compatibility  relation  (i.e.,  Equation  2.37)  along  this  characteristic 
line  can  be  written  as  follows 


(£)■  (p2-p3)  +  u3  =  {-1,3  +  a(T  +  -y)}  4'r  (s-k2) 

ave  '  ave 

where  p  at  ©  is  unknown.  The  line  ©  -  ©  is  a  stream  line  for 
both  gas  and  particle  flows.  Thus  from  Equation  (2.36) 


SQ  -  S.  =  R.  AT 
2  1  4 'ave 


1,3.4  ,) 


and  from  Equation  (2.40) 

Bu 

(pP2  "  PP1)  =  (B2  '  5^  >ave  4T  (£-41,) 

These  equations  along  with  the  definition  of  entropy  determine  p, 
and  S  at  point  ©  .  In  practice,  an  iteration  process  involving  an 
average  speed  of  sound  between  ©  and  ©  is  requii’ed  to  accurately 
locate  the  point  (3)  .  Normally,  this  iteration  converges  in  about 
three  cycles . 
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The  second  stationary  boundary  of  the  problem  is  the  supersonic 
exit  plane  of  the  nozzle  located  at  X  =  1  in  region  [2],  In  single 
phase  flow  (gas  only),  a  common  procedure  is  to  linearly  extrapolate 
all  variables  from  the  two  adjacent  upstream  mesh  points  assuming  they 
are  in  supersonic  flow.  Errors  introduced  by  this  extrapolation  never 
influence  the  upstream  calculations  since  the  flow  velocity  exceeds  the 
local  speed  of  propagation  of  a  disturbance.  This  method  yields  accept¬ 
able  results  in  the  present  analysis  when  the  particle  flow  is  deleted. 
With  the  particle  flow  included,  however,  the  small  errors  in  u^  and 
pp  at  the  boundary  do  propagate  upstream  and  eventually  destroy  the 
solution.  A  speed  of  propagation  characteristic  of  the  suspended 
particles  alone  has  not  been  defined.  Since  the  speed  of  sound  in  a 
gas-particle  mixture  is  less  than  its  value  in  the  gas  alone,  attempts 
to  explain  this  behavior  on  the  basis  of  an  effective  speed  of  pro¬ 
pagation  have  not  been  successful.  The  current  analysis  avoids  this 
difficulty  with  a  complete  method  of  characteristics  solution  at  the 
exit  plane  for  all  five  flow  field  variables.  This  is  illustrated  on 
a  T  -  X  diagram  in  Figure  7  for  the  example  of  Figure  5.  The  compati¬ 
bility  relations  along  these  directions  are  written  in  a  form  similar 
to  Equations  (2.42)  through  (2.44)  and  a  standard  iteration  technique 
is  used  to  locate  the  intersection  points  ©,  0>  ©  and  (o)  .  This 
solution  at  the  supersonic  boundary  eliminates  the  upstream  error  pro¬ 
pagation  observed  with  the  extrapolation  technique. 

The  properties  upstream  and  downstream  of  a  traveling  shock  wave 
are  determined  at  the  next  time  step  in  an  analogous  manner.  For  the 
right-running  shock  wave  in  Figure  5>  upstream  properties  are  located 
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Figure  7.  Method  of  Characteristics  Solution  to 
Obtain  Boundary  Values  at  Supersonic 
Exit  Plane 

at  X  =  0  in  region  [2]  and  downstream  properties  are  at  X  -  1  in  region 
[1].  These  two  points  are  treated  separately  in  the  computational 
method,  but  they  represent  the  same  point  m  physical  -Pace- 


integration  procedure  is  shown  in  the  physical  t  -  x  plane  by  Figure  8. 
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Figure  8.  Method  Used  to  Compute  Shock  Wave  Propagation 
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The  shock  wave  velocity  at  time  t  =  t  is  integrated  with  respect  to 

time  to  give  the  new  location  (same  as  (0)  ).  In  agreement  with 

intuition,  the  method  to  predict  the  upstream  properties  at  this  new 

location  is  based  only  on  information  ahead  of  the  wave.  This  is 

accomplished  with  a  complete  method  of  characteristics  solution  in  the 

triangular  region  ®  -  ©  -  ©  .  The  routine  is  the  same  one  used 

to  obtain  the  supersonic  boundary  condition.  Since  the  suspended  solid 

particles  exert  no  pressure,  they  are  not  influenced  by  the  discontinuous 

jump  through  uhe  shock  wave.  Hence  values  of  p  and  u  determined  in 

P  P 

the  upstream  property  solution  at  are  assigned  to  (£j)  .  Of  course, 

the  particle  flow  is  involved  in  a  relaxation  process  downstream  of 
the  wave.  Three  gas  properties  (p,  u,  S)  at  the  downstream  location  of 
the  shock  wave,  \2y  ,  and  the  new  shock  wave  velocity,  w  ,  remain  to 
be  determined.  Three  of  the  necessary  equations  come  l.-om  the  Rankine- 
Hugoniot  relations  which  conserve  mass,  momentum  and  energy  across  the 
discontinuity.  In  a  stationary  coordinate  system  where  subscripts  1 
and  2  indicate  upstream  and  downstream  respectively,  these  equations 
are 


P2U2  -  Vp2  -  pl>  ‘  P1U1 


(2.45) 


h2^/y  +  Plul^  +  ws) 


PlU2Vs 


=  Pl/Y 


P1U1 


(2.46) 
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(2.47) 


£r  ♦  V«l  -  %) 


The  shock  wave  velocity  w  relative  to  the  fixed  coordinate  system  is 
positive  for  a  right-running  wave.  The  fourth  equation  is  the  compati¬ 
bility  relation  along  the  trailing  subsonic  characteristic  which  reaches 
the  shock  from  the  high  pressure  side.  In  the  present  example,  this  is 
the  right-running  characteristic  -  (^)  .  The  solution  to  this 

system  of  four  simultaneous  equations  is  obtained  with  a  minimization 

(44  45) 

scheme  due  to  Davidon  v  *  .  The  values  of  p,  p  ,  u,  u  and  S  deter- 

p  p 

mined  on  the  downstream  side  of  the  shock  wave  are  assigned  to  the 
variables  at  X  =  1  in  region  [l]  and  the  upstream  values  are  assigned 
at  X  =  0  in  region  [2],  Thus,  all  boundary  conditions  are  obtained  in 
a  manner  consistent  with  the  conservation  equations  and  no  continuous 
differential  equation  is  integrated  across  a  discontinuity. 

It  is  anticipated  that,  under  certain  conditions,  continuous 
disturbances  in  the  rocket  engine  flow  field  could  steepen  and  eventually 
form  shock  waves.  Hence,  use  of  the  above  computational  method  is  con¬ 
tingent  upon  recognizing  the  formation  of  an  embedded  shock  wave  in  a 
region  of  continuous  flow.  Any  attempt  to  equate  an  infinite  pressure 
gradient  (indicating  the  presence  of  a  shock  wave)  with  a  large  rumber 
computed  from  finite-difference  values  is  arbitrary  at  best.  Hence, 
a  more  sensitive  method  is  required  to  predict  the  onset  of  the  dis¬ 
continuity.  Figure  9  illustrates  a  typical  situation  in  the  numerical 
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Figure  9.  Pressure  Waveform  During  Steepening  of  a  Compression 
Wave  in  the  Numerical  Flow  Field 

flow  field  when  a  compression  wave  is  steepening.  Following  Moretti's 
work  ^  ,  a  third-order  polynomial  fitted  to  the  values  at  four  mesh 
points  is  assumed  to  represent  the  pressure  distribution  in  this  interval. 
If  this  polynomial  predicts  an  infinite  gradient,  its  location  is  taken 
as  the  origin  of  the  discontinuity.  Numerically,  this  is  done  in  the 
following  manner.  Given  the  two  functions,  p  =  F(x)  and  x  =  G(p)  in 
the  interval  x^  £  x  £  x^,  the  location  at  which  dF/dx  -*  »  is  also  defined 
by  dG/dp  =  0.  Thus  the  polynomial., 


op 

x  =  G(p)  =  ApJ  +  Bp  +  Cp  +  D 


(2.48) 
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is  specified  by  the  constants  A,  B,  C,  D  determined  by  the  given  values 
at  the  mesh  points,  (p^,  x^)  i  =  1  .  .4.  The  zero  derivative  must 
satisfy 


§=°  =  §=3Ap2  +  2Bp  +  C 


(2.49) 


Real  roots  of  Equation  (2.4-9)  are  possible  only  for  non-negative  values 
* 

of  the  discriminant  D,  where 


(2.50) 


In  the  present  problem,  D  maintains  a  negative  value  for  continuous 
pressure  distributions  which  do  not  steepen.  This  provides  an  easy 
check  for  the  steepening  of  compression  waves.  At  each  time  step,  a 
computer  search  routine  checks  for  a  non-negative  value  of  D  for  all 
combinations  of  four  adjacent  mesh  points  in  the  continuous  regions; 
if  one  is  found,  the  x  location  of  the  roots  of  Equation  (2.49)  are 
calculated  with  Equation  (2.48).  If  a  root  lies  within  the  interval 
x^  s  x  s  x^,  an  interior  boundary  is  inserted  at  this  location  and 
the  computation  proceeds  tracking  i  Mach  wave  as  a  discontinuity. 
Hence,  as  the  compression  wave  continues  to  steepen,  it  is  treated 
as  a  shock  wave. 

It  might  be  argued  that  this  "sensing"  method  is  heavily 
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dependent  on  how  well  the  third-order  polynomial  represents  the  actual 
pressure  distribution.  In  spite  of  its  potential  short-comings,  the 
technique  has  proved  quite  sensitive  to  discontinuity  formation  in  the 
present  investigation.  Even  if  the  formation  is  predicted  prematurely, 
there  are  no  penalties  associated  with  tracking  a  Mach  wave  as  a  dis¬ 
continuity.  On  the  other  hand,  the  author  can  verify  that  severe 
penalties  are  incurred  for  late  recognition  of  discontinuous  behavior. 

Integration  of  the  equations  in  the  floating  computational 
coordinate  system  with  moving  internal  boundaries  requires  many  decisions 
of  a  practical  nature.  Typical  of  these  is  how  a  shock  wave  passes 
through  the  exit  plane  of  the  nozzle.  As  the  right-running  shock  wave 
in  the  example  of  Figure  5  travels  through  the  throat  and  toward  the 
exit  plane,  the  physical  length  c  -  b  of  region  [2]  becomes  small.  An 
infinite  number  of  calculations  would  be  necessary  to  reach  the  exit 
plane  since  the  forward  marching  time  step  (see  next  section)  is  pro¬ 
portional  to  this  length.  In  the  present  analysis,  when  c  -  b  of 
region  [2]  is  less  than  one  percent  of  the  total  flow  field  length, 
region  [2]  is  deleted  from  the  computation.  The  downstream  boundary 
of  region  [l]  is  adjusted  to  the  exit  plane;  hence,  the  shock  wave  has 
left  the  system.  Repeating  this  calculation  with  a  criterion  of  0.1 
percent  has  no  effect  on  the  results  other  than  to  lengthen  the  com¬ 
putation  time. 


Forward  Marching  Stability 

Explicit  integration  techniques  applied  to  partial  differential 
equations  must  observe  a  forward  marching  stability  criterion.  Other- 
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wise,  small  errors  associated  with  the  finite  difference  approximation 
grow  unbounded  as  the  solution  progresses.  The  well-known  von  Neumann 
stability  analysis  uses  the  linearized  equations  topredict  the 
growth  or  decay  of  these  error  waves.  However,  the  analysis  is  impracti¬ 
cal  to  perform  for  the  present  system  of  equations  and  the  result  would 
be  strictly  valid  only  for  the  linearized  form.  The  success  of  many 

numerical  solutions  employing  the  Courant-Friedrichs-Lewy  (C.F.L. ) 

(47) 

Condition  indicates  that  it  is  an  acceptable  substitute.  The 
C.F.L.  Condition  requires  the  domain  of  dependence  of  each  point  in 
the  finite-difference  flow  field  to  include  its  physical  domain  of 
dependence.  This  restriction  on  the  step  size  means  that  the  speed  of 
propagation  of  numerical  information  will  equal  or  exceed  the  speed  of 
propagation  in  the  physical  flow.  This  can  be  seen  with  the  aid  of 
Figure  10  where  the  point  B'  is  located  at  the  intersection  of  the 
characteristic  line  with  the  maximum  slope  through  point  A  and  its 
mirror  image  through  point  C.  With  the  maximum  slope  defined  as  |u| 

+  a,  the  C.F.L.  Condition  requires 


Ax 

At 


£  u  +  a 


(2.51) 


or 


At  £ 


Ax 

u  +  a 


(2.52) 
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Figure  10.  Illustration  of  the  C.  F.  L.  Condition 


The  present  analysis  uses  this  stability  condition  in  the  form. 


At 


minimum  over 
the  flow  field 


(2.53) 


where  the  value  of  the  constant  was  influenced  by  Rentier's  ^  ^  analysis 
of  group  velocity  and  propagation  of  numerical  errors. 

It  is  sometimes  reported  that  a  well-behaved  solution  to  a  non¬ 
linear  problem  was  gradually  or  even  suddenly  destroyed  by  numerical 
wiggles.  Their  appearance  is  often  attributed  to  "non-linear  instability". 
The  present  investigation  subscribes  to  Moretti's  theory  that 

numerical  non-linear  instability  does  not  exist.  Deterioration  of  a 
numerical  solution  can  always  be  traced  to  an  improper  numerical  treat- 
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ment.  Computations  which  integrate  continuous  equations  across  dis¬ 
continuous  phenomena,  incorrectly  determine  boundary  values,  or  fail 
to  recognize  regions  of  large  truncation  error  will  generate  non¬ 
physical  solutions  to  the  equations. 

Recently,  Cheng  '  has  questioned  the  ability  of  an  explicit 
integration  technique  to  accurately  predict  oscillatory  behavior.  Part 
of  the  discussion  is  directed  toward  a  monotonically  damped  wave  solution 
of  Burger's  Equation,  i.e., 


u,  +  uu  -  u  =0 

t  X  XX 


Using  the  two-step  Lax-Wendroff  technique  (equivalent  to  a  second- 

order  Taylor  series)  and  evaluating  initial  conditions  and  boundary 
values  from  a  known  analytical  solution,  Reference  (36)  shows  a  pro¬ 
gressively  deteriorating  approximation  to  the  known  solution.  The 
solution  is  no  longer  useful  after  a  computation  of  one  wave  length. 

The  failure  is  attributed  to  phase  errors  caused  by  dispersion,  with 
the  suggested  cure  being  difference  schemes  of  higher  than  second-order 
accuracy.  This  has  a  serious  impact  on  the  present  combustion  instability 
analysis  which  is  required  to  predict  oscillatory  motion,  with  or  with¬ 
out  shock  waves. 

A  simple  test  case  was  constructed  from  the  general  equations 
to  examine  this  difficulty  in  detail.  For  small  amplitude  disturbances, 
the  dependent  variables  can  be  expressed  in  the  form, 
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p(x,t)  =  p(x)  +  p’(x,t) 

p(x,t)  =  p(x)  +  p'(x,t) 

u(x,t)  =  U(x)  +  u*(x,t) 

where 

(p'/p)2  «  1  etc. 

When  combustion,  solid  particles  and  mean  flow  are  absent,  Equations 
(2.9)  and  (2.1l)  reduce  to, 

Continuity 
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2 

where  p'  =  p'/ya  follows  from  isentropic  flow.  This  system  is  equi¬ 
valent  to  the  linear  second  order  wave  equation.  The  behavior  of  simple 
wave  motion  between  two  solid  boundaries  was  computed  as  a  function  of 
time.  Following  Reference  (36),  an  exact  analytical  solution  was  used 
to  evaluate  the  initial  conditions  and  all  boundary  values .  After  three 
cycles,  the  numerical  solution  no  longer  resembled  the  initial  wave  form. 
However,  re-running  the  same  problem  using  the  method  of  characteristics 
to  determine  boundary  values  resulted  in  a  wave  form  that  was  indis¬ 
tinguishable  from  the  initial  wave  after  23  cycles.  This  behavior  can 

be  explained  by  the  fact  that  the  analytical  solution  to  a  set  of 

(!<g) 

differential  equations  is  not  the  same  ^  as  the  analytical  solution 
to  the  corresponding  finite-difference  equations.  If  the  equations 
describe  oscillatory  motion,  there  is  a  small  phase  difference  between 
the  corresponding  solutions  which  grows  in  time.  Thus,  using  the  exact 
analytical  solution  of  the  differential  equations  to  specify  boundary 
conditions  to  a  numerical  solution  of  the  finite-difference  equations 
means  that  a  slight  error  is  committed  at  the  boundary.  These  errors 
propagate  inward  and  eventually  destroy  the  solution. 

It  is  concluded  that  forward  marching  integration  based  on  a 
second-order  Taylor  series  expansion  which  observes  the  Courant -Friedrichs - 
Lewy  Condition  and  obtains  boundary  values  with  the  method  of  character¬ 
istics  is  internally  consistent  and 'hence  stable. 
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CHAPTER  III 


THE  COMBUSTION  RESPONSE  MODEL 

Problem  Formulation 

The  response  model  of  the  burning  solid  propellant  must  predict 
how  the  combustion  process  is  altered  by  changes  in  pressure  in  the  thin 
flame  zone  adjacent  to  the  ,  urface.  This  information  is  supplied  to  the 
combustion  chamber  analysis  of  Chapter  II  through  the  time -dependent 
quantities  m^,  m^,  and  hp  which  enter  the  control  volume  perimeter 
boundary  from  the  propellant  flame  zone.  The  combustion  process  which 
specifies  this  variation  is  strongly  dependent  on  the  nature  of  the 
disturbance  in  the  chamber,  suggesting  the  necessity  of  a  simultaneous 
solution  to  both  problems.  This  will  be  done  in  the  present  study  in 
the  following  manner. 


-x 

■y 


Figure  11.  Method  Used  to  Represent  Cylindrically 
Perforated  Solid  Propellant 


73 


It  will  be  assumed  that  the  combustion  of  the  cylindrically  perforated 
propellant  can  be  represented  by  a  series  of  fixed  burning  stations  as 
shown  in  Figure  11.  In  the  flame  zone  adjacent  to  the  propellant, 
property  variations  in  the  y  direction  are  much  greater  than  those  in 
the  x  direction,  which  allows  the  process  to  be  treated  as  one-dimensional 
at  each  station.  Each  local  burning  solution  will  depend  in  the  chamber 
flow  conditions  in  the  control  volume  directly  below,  i.e.,  at  the  same 
axial  location.  In  this  way,  the  combustion  response  of  the  solid 
propellant  to  arbitrary  disturbances  in  the  chamber  is  determined  by 
the  history  stored  at  each  location,  over  the  length  of  the  grain. 

The  combustion  model  used  at  each  burning  station  along  the 

propellant  is  based  on  four  simplifying  assumptions  which  have  been 

(26) 

employed  in  nearly  all  combustion  instability  studies  to  date  '  . 

These  assumptions  are: 

(1)  The  solid  phase  (unburned  propellant)  is  homogeneous  and 
one -dimensional.  Thus,  the  analysis  will  be  unable  to 
distinguish  between  homogeneous,  double-based,  and  composite 
propellants.  The  solution  must  be  based  on  average  values 
of  propellant  composition.  Three-dimensional  surface 
topological  effects  are  lost  in  the  one-dimensional  approxi¬ 
mation  which  means  that  the  surface  regression  rate  must 

be  viewed  as  an  an  average  rate. 

(2)  The  conversion  of  the  solid  phase  to  gas  is  represented 
by  a  simple  pyrolysis  law.  Usually  the  interface  region 
between  solid  phase  and  the  gas  phase  is  collapsed  to  a 
plane,  which  effectively  "lumps"  all  decomposition, 
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pyrolysis  and  other  surface  reactions  into  a  single  reaction. 
This  overall  pyrolysis  reaction  is  taken  to  be  an  Arrhenius 
type,  which  may  have  either  an  endothermic  or  exothermic  net 
heat  release.  The  present  analysis  will  assume 


*  * 
r  =  regression  rate  =  Br  e 


* 

T 

o  sur 


where  the  possible  pressure  dependence  of  B*  is 

neglected. 

(3)  Condensed  phase  reactions  are  neglected.  In  view  of 

Prices  description  of  the  combustion  of  metal  loaded 

solid  propellants,  this  assumption  necessarily  obscures 
important  details.  Except  for  the  exp3.oratory  work  of 
Kumar  and  Culick  with  a  global  condensed  phase 
reaction,  this  assumption  is  state-of-the-art.  Without 

a  proper  description  of  the  generation  of  solid  oxide 

particles,  the  present  analysis  will  assume  that  the  mass 

flow  rate  of  particles  to  the  control  volume,  ,  is  a 

constant  fraction  of  the  gaseous  mass  flow  rate,  m  .  This 

6 

fraction  is  taken  as  the  percentage  of  metal  loading  in 
the  unburned  propellant. 

(4)  The  behavior  of  the  gas-phase  flame  zone  is  quasi-steady. 
The  specific  heats  of  the  gas-phase  constituents  are  the 
same  order  of  magnitude  as  the  specific  heat  of  the  solid 
propellant.  Since  the  density  of  the  unburned  solid 
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propellant  Is  very  much  greater  than  the  density  of  the  gas- 
phase,  the  volumetric  heat  capacity  of  the  solid  compared 
to  the  gas  is  of  the  order  of  the  density  ratio.  This  fact 
supports  the  argument  that  the  gas-phase  will  adjust  to 


thermal  changes  much  more  rapidly  than  the  solid  phase, 


i-e.,  ■& 


«xrl. 


&t 'gas-phase  bt  solid-phase’ 


In  addition,  the 


characteristic  times  for  most  chemical  reactions  are  much 


smaller  than  thermal  wave  propagation  times  in  the  solid. 
Thus,  it  is  assumed  that  the  controlling  time -dependent 
process  is  the  thermal  wave  propagation  in  the  solid,  and 
the  flame  region  can  be  described  by  a  steady-state  solution 
adjusted  at  each  instant  of  time.  This  affords  a  great 
simplification  if  the  flame  temperature  is  assumed  known, 
since  the  analysis  becomes  insensitive  to  the  time-dependent 
kinetics  and  transport  processes  in  the  gas-phase  flame. 

The  accuracy  of  the  quasi-steady  assumption  is  still 
open  to  question.  It  is  often  claimed  that  this  approxi¬ 
mation  should  remain  valid  for  cases  of  low  or  intermediate 


frequency  instabilities  such  as  found  in  the  present  study. 
However,  a  detailed  and  sophisticated  analysis  will  be 
required  to  verify  this .  Further  comment  will  be  made  in 
Chapter  V. 

Mary  well-known  combustion  response  theories  based  on  pre;  sure 

(nC  \ 

coupling  have  been  compared  in  an  excellent  review  by  Culick  .  All 
contain  the  standard  assumptions  outlined  above,  although  they  differ 
in  the  treatment  of  details.  Each  theory  has  assumed  the  existance  of 
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small  amplitude  periodic  oscillations  in  the  gas-phase  flame  zone.  It 
is  shown  that  as  a  consequence  of  the  above  assumptions,  each  linearized 
analysis  predicts  the  same  two-parameter  (A,  B)  form  of  the  combustion 
response  function,  i.e., 


ft  H 


m'/m 

P'/P 


_ nAB _ 

X  +  X/A  -(1  +  A)  +  AB 


(3.1) 


where  m*  =  surface  mass  flow  rate  perturbation 

p'  =  pressure  perturbation 

n  =  pressure  index  for  time-independent  burning 

X  =  function  of  frequency  of  oscillation  (complex  number) 

The  steps  leading  to  Equation  (3.1)  can  be  found  in  Reference  (26). 

This  equation  predicts  the  trend  of  experimental  data  for  acoustic 

(27) 

oscillations  in  a  T-burner  ,  etc.  A  plot  of  ft  versus  fl,  defined 
as  the  ratio  of  frequency  to  burning  rate  squared,  shows  evidence  of 
a  resonance  effect  at  a  non-zero  value  of  frequency,  determined  by 
the  values  of  A  and  B  (see  Figure  12).  The  shape  of  the  curve  and 
the  peak  response  varies  with  the  values  of  A  and  B.  However,  even 
with  the  adjustable  parameters,  Equation  (3.1)  ij  not  able  to  correlate 
all  available  experimental  data.  The  degree  of  simplicity  which  can 
be  tolerated  in  a  theoretical  model  of  the  combustion  process  remains 
an  open  question. 

(34) 

The  instability  analysis  of  Levine  and  Culick  uses  Equation 
(3.l)  as  the  basis  for  the  transient  combustion  response.  Since  ft  is 
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Figure  12.  Typical  Two -Parameter  Combustion  Response 
Curves  as  a  Function  of  Frequency 

a  given  function  of  frequency,  an  inverse  Laplace  transform  must  be 
evaluated  to  transform  the  results  into  the  real-time  plane.  It  is 
not  clear,  however,  how  results  derived  for  small  amplitude  periodic 
oscillations  can  be  applied  in  an  analysis  of  the  transient  combustion 
response  which  may  accompany  arbitrary  amplitude  disturbances.  To 
avoid  this  difficulty,  the  present  analysis  will  solve  the  full  non¬ 
linear  time -dependent  equations  of  the  combustion  model  simultaneously 
with  the  conservation  equations  in  the  combustion  chamber. 

The  combustion  response  model > in  the  current  study  diverges 
from  previous  efforts  in  two  areas;  these  are: 

(l)  It  will  not  be  assumed  that  a  disturbance  in  the  gas -phase 
is  oscillating  with  a  known  frequency,  and  instead,  the  full  transient 
behavior  of  the  combustion  response  will  be  determined. 
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(2)  No  restrictions  will  be  placed  on  the  magnitudes  of  the 
gas  and  solid  phase  disturbances. 

The  mathematical  model  of  the  solid  propellant  combustion  process 
follows  from  a  development  in  Reference  (3*0  and  employs  the  four 
assumptions  outlined  in  the  beginning  of  this  chapter.  In  Figure  13 
is  a  description  of  the  problem  as  seen  from  a  coordinate  system 
attached  to  the  propellant  surface,  which  is  regressing  at  speed  r*. 


Figure  13.  One -Dimensional  Model  of  the  Solid 
Propellant  Combustion  Process 

At  the  chamber  pressures  involved  in  the  present  study,  both  diffusion 

and  chemical  kinetics  will  influence  the  structure  of  the  flame  zone. 

Since  a  detailed  gas-phase  solution  is  beyond  the  scope  of  this  work, 

the  combustion  process  is  assumed  to  be  distributed  in  such  a  manner 

that  the  flame  zone  exists  from  y*  =  0  to  y*  =  y*.  This  is  in  contrast 

to  a  flame  sheet  approximation  (i.e.,  all  combustion  occurs  at  y*  =  y*^,) 

which  may  be  valid  at  low  pressures  where  reaction  rates  are  much 
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slower  than  mixing  processes.  A  simple  one-step  global  reaction  is 
assumed  in  the  flame;  i.e., 

Oxidizer  +  Fuel  -*  Products 

*  *2 

implying  a  gas-phase  reaction  rate,  cu  ,  proportional  to  p  .  It  is 
further  assumed  that  the  product  of  reaction  rate  (actually,  production 

-ft 

rate  of  products)  and  heat  release,  cu  Q^,  can  be  taken  as  constant 
*  “ 

across  the  width  of  the  flame.  This  set  of  assumptions  permits  a 

(33) 

separation  of  the  problem  into  three  distinct  regions  '  ,  each  of 

which  is  analyzed  individually.  These  are: 

* 

(l)  The  non-reacting  solid  phase  (-  <=  £  y  £  0)  region.  The 
unburned  solid  propellant  is  assumed  to  be  incompressible  and  hence  the 
equations  of  continuity  and  momentum  become  trivial.  The  energy 
equation  can  be  written  as. 


*  *  dT*  .  *  *  *  ST*  *  d2T*  . 

p  c  — s  +  p  c  r  — ;  -  k  — =  0 

Ks  s  *  Fs  s  ,*  s  *2 

ot  oy  oy 


(3.2) 


If  the  overall  pyrolysis  of  solid  to  gas  follows  a  first  order  Arrhenius 
type  reaction,  then 


* 


r 


* ,  * 
„  -E  /R  T 

*  S3'  n  ! 

B  e 
r 


0  sur 


(3.3) 
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where 


* 

B 

r 


=  constant 


T  =  surface  temperature 
sur 

(2)  The  solid-gas  interface  (y  =  0).  For  a  control  volume 

* 

with  zero  thickness  which  surrounds  the  plane,  y  =  0,  the  conservation 
of  mass  yields, 


m  =  surface  mass  flow  rate  = 


*  * 
P  r 


o 


*  *. 
pgvgl© 


(3.4) 


A  first  integral  of  the  energy  equation  evaluated  between  the  pro¬ 
pellant  side  (-)  and  the  gas  flame  side  (+)  gives 


*  * 
m  h 

s 


0 


*  * 
—  th  h 

g 


© 


-D 


*  st 

V  * 
s  3y 


* 


(3.5) 


Letting  Q  be  the  heat  of  reaction  for  the  surface  pyrolysis  reaction, 
s 


© 


(3.6) 
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where  Q,  is  positive  for  an  endothermic  heat  release,  Equation  (3. 5) 
s 

can  be  rewritten  as, 


L  s  Sy  Je 


*  * 

=  -  m  Qs  + 


K  5  ] 


by  © 


(3.7) 


This  relationship  is  the  basis  for  the  time -dependent  boundary  con¬ 
dition  imposed  on  the  thermal  profile  in  the  unburned  solid. 

(3)  The  gas-phase  flame  zone  (0  £  y  £  yf).  When  the  assumptions 
of  a  thin  flame  zone  and  quasi-static  behavior  (i.e.,  -*  o)  are  applied 
to  the  conservation  equations  for  mass  and  momentum,  the  result  is. 


*  *  .* 
p  v  =  constant  =  m 
g  g 


(3.8) 


* 

p  =  constant 


(3-9) 


These  constant  values  across  the  flame  zone  are  the  instantaneous 
values  of  mass  flow  rate  and  pressure,  and  hence  they  vary  with  time. 
The  quasi-steady  energy  equation  can.be  written  as 


.*  *  bT 

m  — * 

P  -s  * 

Sy 


9  * 

*  *  * 
kg  *2  “  Qf  ® 

oy 


(3.10) 
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where 


*  * 

constant  0  £  y  s 

-X-  -x* 

0  yf  <  y  <  +  ro 


Again,  since  the  reaction  rate  is  a  function  of  the  pressure  in  the 


flame,  the  "constant"  value  of  Q„cu  changes  with  time. 


In  terms  of  the  following  non-dimensional  variables, 


*  * 

Qfco  = 


*  -  * 

T  =  T  /T 

*  ,  ** 
y  s  y  Ayc 

t  =  t  (a r^  ) 


I 

i 


*  *A 
r  =  r  /a  y 

/  i 


r“  c 


*  *  *«  * 
a  s  *  cpL  yc/kg 


(3.11) 


* .  *  * 
z,  =  ?  /c  T 
1  s'  s  r 
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where 


cu  (t)  = 


*.  * 
2  Wf 
“Wp  e 


and 


the  equations  in  the  three  regions  can  be  rewritten  as  follows 
(l)  Solid  phase  (-  <*>  <  y  s  0) 


T,  +  rT 
t  y 


T  =  0 

yy 


(2)  Interface  (y  =  0) 


T 


y 


r 


(3.12) 


(3.13) 
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(3)  Gas  phase  (0  £  y  <  +  ») 


where 


°'Ty  '  Tyy  =  5  <3.14) 


A 


Q  = 


constant  0  s  y  5  y^, 


o  yf  <  y  <  +  ro 


The  temperature  distribution  in  the  flame  zone  can  be  found 
immediately  if  it  is  assumed  that  the  gases  achieve  the  flame  temperature 
when  combustion  is  complete.  Thus,  if  T  =  T  f  as  y  -♦  +  <»,  Equation 
(3.14)  yields 


1  (y)  =  T&f  for  (yf  s  y  <  +  »)  (3-15) 


Matching  the  .emperature  and  its  first  derivative  at  y  =  y^.,  the  thermal 
profile  in  the  fiime  zone,  0  £  y  £  y^.,  must  be 


A  A  Qf(y-yf.)-1 

T(y)  =  Taf  +  ^  (y  -  yf)  +  4  [1  "  e  ] 


(3.16) 


O' 
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(0  S  y  S  yf) 

The  temperature  gradient  in  the  gas  phase  at  the  interface  is  then  given 

by, 

TyUo  “  ^  Cl  '  e'Wf  ]  (3'17) 

Since, 

*  *  *  *  * 
cu  yf  =  m  =  pgr 

it  can  easily  be. shown  for  the  cases  of  interest  in  the  present  study 
that 

2 

oyf  «  lOr 

This  means  that  the  value  of  the  exponential  term  in  Equation  (3.17) 
is  negligible  for  all  but  very  low  burning  rates  (see  also  comments 
in  Reference  3*0-  Hence,  the  temperature  gradient  of  the  gases  at 
the  interface  may  be  taken  as, 
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V©  =TyU  =  ,/“ 


(3.18) 


Equation  (3.13)  can  then  be  written  as, 


where 


ryl  0  ’  ’  V  + 


(3.19) 


%  “  z2(t)  =  (  h  )  (  -TT«  )  (  4-»  ) 

'  ks  ;  c  T  7  p  a  7 
p  r 


p  a 
Ks  r 


This  equation  is  the  time -dependent  heat  transfer  boundary  condition 
imposed  on  the  temperature  profile  in  the  solid.  Thus,  the  system 
reduces  to  the  problem  of  describing  the  behavior  of  the  temperature 
distribution  in  the  unburned  solid  propellant,  which  is  governed  by 
the  following  init ia.L -value  problem: 


Tt  +  rTy  -  T^  =  0  (-sys  0)  (3.20) 


T(y,0)  =  To(y) 


(a) 


I.C. 
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B.C. 


(b) 


Ty(0,t)  =  -Z^r  +  Z2(Qp/r  =  6(t)  (c) 


due  to  due  to 
reaction  reaction 
at  surface  in  flame 


The  instantaneous  value  of  the  surface  temperature,  Tgur  =  T(o,t)  , 
predicted  by  this  solution  determines  the  surface  regression  rate  r, 


r 


=  B  e 
r 


-e*/R  t* 

s'  o  sur 


(3-21) 


which  specifies  the  mass  flow  rate  from  the  surface  according  to 


m  -  <psVVr 


(3.22) 


Since  this  theory  assumes  a  homogeneous  solid  phase,  it  is  unable  to 
distinguish  between  mass  flow  rates  of  gas  and  particles.  The  present 
analysis  assumes  that  the  weight  percent  leading  of  solids  in  the 
unburned  propellant  can  be  used  directly  to  determine  m^  and  m^  which 
enter  the  control  volume  perimeter. 

Method  of  Solution 

The  solution  of  Equation  (3.20)  is  obtained  with  a  method  of 


88 


invariant  imbedding  discussed  by  Meyer  where  the  desired  solution 

is  imbedded  in  a  family  of  initial  value  solutions  as  opposed  to  a 

(52) 

family  of  boundary  value  solutions  (e.g.,  see  Lee'  ') .  Meyer's 

(53) 

method  follows  from  the  early  work  of  Vichnevetsky  .  This  approach 
makes  use  of  the  "method  of  lines,"  whereby  a  partial  differential  equation 
is  converted  to  an  ordinary  differential  equation  by  discretizing  all 
but  one  of  the  independent  variables,  which  is  left  continuous.  Thus, 
Equation  (3 .20)  becomes  an  ordinary  differential  equation  in  the  space 
variable,  y,  after  discretizing  the  time  derivative  with  a  backward 
difference.  Treating  this  equation  as  an  initial-value  problem,  the 
general  solution  can  be  found  with  the  method- of- characteristics .  Imposing 
the  boundary  conditions  determines  the  desired  solution  at  the  current 
time  level.  Since  this  technique  is  impl.‘ cii ,  it  is  stable  in  a  forward 
marching  time  integration  and  is  not  sat  ,eci  to  the  usual  stepsize  stability 
control  required  in  an  explicit  calcui  .tion.  Of  course,  excessive  values 
of  the  incremental  time  step  lead  to  unacceptable  truncation  error. 

For  convenience,  use  the  notation. 


$  =  l(y)  H  T(y,  time  "n") 


=  l,n"1(y)  =  T(y,  time  "n-1") 


where  all  information  is  assumed  known  at  time  "n-1."  Discretize  the 
partial  time  derivative  in  Equation  (3 .20)  with  a  backward  difference  to 


obtain, 
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(3.23) 


Tt  =  (T  -  Tn_1)/At 


Substituting  Equation  (3-23)  into  Equation  (3-20)  results  in  a  second, 
order  ordinary  differential  equation  in  y  at  time  V\  viz. 


T 


"  =  rT'  +  At 


T* 


1_  An-1 
At 


(3.24) 


where 


The  boundary  conditions  to  Equation  (3.24)  follow  directly  from  Equations 
(3.20(b))  and  (c), 


T(-y0) 


(3.25) 


T'(o)  =  G  =  G(time  "n") 


(3.26) 


To  circumvent  a  numerical  integration  to  minus  infinity,  the  solid 
boundary  condition  is  applied  at  the  finite  location,  -  y0-  Obviously 
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this  point  must  be  located  deeper  in  the  solid  than  any  thermal  wave 
can  penetrate. 

By  defining  the  variable,  v,  as 


v 


(3.or?) 


Equation  (3*24)  can  be  represented  by  the  equivalent  first  order  system, 


(3.28) 


v' 


rv  + 


JL_ 

At 


T 


_1_  £n-l 
At 


(3-29) 


If  these  equations  are  regarded  as  characteristic  equations  for  an 

(54) 

initial  value  problem  then  characteristic  theory  for  a  first  order 
partial  differential  equation  states  that  the  solution  of  Equations 

A 

(3.28)  and  (3.29)  generates  the  integral  surface  T(y,v)  of  the  equation 
(see  Appendix  A), 


ST  or  r 

57 +  ~ Crv 


ST 

Sv 


—  T  -  —  Tn_1]  = 
At  At  J 


(3.30) 
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The  general  solution  of  Equation  (3-30)  can  he  shown 
form, 


to  have  the 


(51) 


T(y,v)  =  e(y)v  +  w(y) 


Substitution  of  Equation  (3.3l)  into  (3-30)  shows  that  0(y) 
must  satisfy  the  equations 


e*  =  i 


re  -  sr  «2 


v.  =  -  i  e(w  -  i"'1) 


and  the  boundary  conditions, 


e(-y0)  =  o 


w(-y  )  =  T 
'  Jo/  cs 


(3-31) 

and  w(y) 

(3.32) 

(3-33) 

(3.34) 

(3.35) 


which  satisfy  Equation  (3.25).  The  solutions  of  Equations  (3-32)  and 

A 

(3*33)  specify  the  integral  surface  T(y,v)  in  which  the  desired  solution 


is  imbedded.  This  particular  solution  is  determined  by  satisfying  the 


remaining  boundary  condition  at  y  =  0;  i.e., 


v|0=G 


(3.36) 


The  solutions  0(y)  and  w(y)  can  be  found  without  knowledge  of  v  by 
integrating  Equations  (3.32)  and  (3-33)  from  (-y  )  to  (0).  Then  the 
above  boundary  condition  (Equation  3-36)  must  be  satisfied  at  y  =  0. 

A 

The  regression  rate  r,  and  G  which  is  a  function  of  r,  depend  on  the 
surface  temperature 


T(0,v|0)  =  6(0)  v|Q  +  w(0)  (3.37) 

This  means  an  iteration  procedure  is  required  to  satisfy  Equations 
(3.32)  and  (3.33)  simultaneously  with  the  boundary  condition,  Equation 
(3.36).  Once  the  boundary  condition  has  been  satisfied,  the  functional 
form  of  v  is  found  from  Equation  (3.29)  integrated  from  (0)  -*  (-yQ), 


v’  -  v[r  +  0(y)/At]  +  [w(y)  -  Tn"1]  (3-38) 


Finally,  the  temperature  profile  is  given  by 
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(3.39) 


T(y,  time  "n")  =  T(y)  =  S(y)v(y)  +  w(y) 


The  advantage  offered  by  this  method  is  the  minimum  amount  of 
computation  time  required.  For  the  special  case  when  r  in  Equation 
(3.24)  is  constant,  only  two  integration  passes  are  needed  to  completely 
determine  the  solution,  i.e.,  (l)  Equations  (3-32)  and  (3*33)  are 
integrated  from  (-y  )  to  (0),  and  the  boundary  condition  Equation  (3 .36) 
is  satisfied  at  y  =  0,  (2)  Equation  (3-38)  is  integrated  from  (o)  to 
(-yQ).  When  the  r  in  Equation  (3.24)  is  not  constant,  a  simple  but 
rapidly  converging  iteration  cycle  is  required  as  mentioned  above.  The 
alternative  way  to  solve  this  second-order  differential  equation  with 
split  boundary  conditions  is  known  as  the  "shooting"  technique.  Begin¬ 
ning  at  one  end  of  the  domain,  the  unknown  value  of  the  dependent 
variable  or  its  derivative  must  be  estimated,  which  allows  the  equation 
to  be  integrated  to  the  opposite  end.  If  the  boundary  condition  at 
this  end  is  not  satisfied,  the  estimated  boundary  condition  must  be 
revised  and  the  integration  repeated.  When  the  equation  or  the  boundary 
conditions  are  nonlinear,  this  process  is  at  best  uncertain  and  may 
require  a  large  number  of  iterations . 

In  the  computer  program,  the  solutions  to  Equations  (3.33)  and 
(3.38)  are  obtained  with  a  numerical  integration  scheme  based  on  the 
trapezoidal  rule.  Equation  (3-32)  is  a  Riccati  equation  which  with 
the  boundary  condition  Equation  (3-34)  can  be  solved  analytically  to 


give, 
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(3.40) 


where 


0(y)  =  [1  -  e 


h^oK/r,  Jd(y+yo) 


3/[6l 


-  ^ 


(-y0  *  y  *  °) 


«i ■  V  ( § )  +  s 


h'-Uiif-k'rt 


?d  ~  ^2  “ 


The  above  algorithm  can  be  summarized  as  follows: 

(1)  Tn_1(y)  is  known;  let  r  =  r[Tn  1(o)] 

(2)  Evaluate  0(y)  from  Equation  (3-40) 

(3)  Solve  the  initial  value  problem 


-  -  ZT  e[w  -  ?'1]  ;  v(-  yo>  =  Tcs 


from  y  =  -y  to  y  =  0 

(4)  Determine, 
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t(o)  =  e(o)  g  +  w(o) 


A 

If  r[T(0)]  differs  from  previous  r,  return  to  (2);  otherwise, 
proceed  to  (5) 

(5)  Integrate  the  equation. 


v'  =  v(r  +  0/At)  +  (w  -  Tn_1) 


from  y  =  0  to  y  =  -  yQ,  using  the  initial  condition  v(o)  = 

A 

G 

(6 )  Evaluate 


T(y)  =  9(y)v(y)  +  w(y) 


(7)  Set 


Tn_1(y)  =  T(y) 


and  return  to  step  (l) 
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Comparison  with  Previous  Models 

With  the  development  of  the  combustion  model  complete,  it  is 

possible  to  compare  the  final  formulation  in  the  present  'study  with  the 

solid  propellant  combustion  analyses  of  the  Princeton  group,  led  by 

(55-59) 

Summerfield.  These  investigations  '  are  directed  toward  ignition 

and  depressurization  problems  in  solid  propellant  rockets,  but  the 
combustion  process  is  treated  in  a  similar  manner.  Each  of  these 
studies  incorporates  the  standard  assumptions  (l)  through  (4)  (see 
pages  56-58)  o.r  their  equivalent.  Hence  the  description  of  the  com¬ 
bustion  response  reduces  to  a  thermal  theory  (as  in  the  present  investi¬ 
gation)  which  determines  a  time -dependent  temperature  distribution  in 
the  unburned  solid.  The  basis  for  comparison  will  be  the  type  of 
initial  value -boundary  value  problem  which  is  derived  in  the  solid 
phase  region,  i.e.,  the  equivalent  of  Equation  (3.20).  Since  the 
analyses  differ  as  to  coordinate  system,  nomenclature,  and  non-d5:.en- 

sionalization,  the  comparison  in  terms  of  symbols  will  be  approximate. 

(55) 

The  work  of  Most,  et.  al.  is  an  ignition  study  designed  to 

predict  the  "start-up"  thrust  transient  in  a  solid  propellant  rocket 
engine  with  low  internal  gas  velocities  (e.g.,  large  port  to  throat 
area  ratio).  Spatial  variations  of  static  pressure  and  temperature 
in  the  chamber  are  neglected.  Ignition  of  the  propellant  is  assumed 
to  occur  when  the  surface  temperature  exceeds  a  given  constant  value. 
Above  this  temperature,  the  burning  rate  is  given  by  the  steady  state 
expression,  r  =  apn.  Below  this  temperature,  the  solid  undergoes  a 
heating  process  due  to  the  ignitor  gases.  The  heating  process  is 


governed  by, 
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T .  -  T  =0 

t  yy 


i.c. 


T(y,0)  =  To(y) 


B.C.  T(-«,t)  =  T 

v  ’  '  cs 

kT  (0,t )  «  q  +  h  (T  -T) 

y  surface  con  gas  sur 

J"  "preignition"  T  ["convective  heat  transfer 
Lsurface  reactions_TLfrom  ignitor  gas 


Since  the  present  analysis  is  concerned  with  the  transient  behavior  in 
the  engine  after  the  propellant  has  ignited,  this  heating  effect  was 
not  studied. 

Krier,  et.  al.  were  the  first  to  present  nonlinear  results 

for  the  unsteady  burning  of  solid  propellants.  Their  analysis  makes 
use  of  the  dependence  of  the  steady  state  regression  rate  on  pressure 
as , 


r  = 


and  on  temperature  as 
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r  =  b(T 


BUT 


T  )m 
cs' 


For  values  of  m  on  the  order  of  10,  the  latter  expression  approximates 
the  temperature  dependence  of  an  Arrhenius  type  law.  Using  these 
steady  state  relations  to  estimate  instantaneous  burning  rates,  the 
time  dependent  problem  for  the  temperature  field  in  the  unburned  solid 
is  given  by, 


T.  +  rT 
t  y 


-  T  = 

yy 


i.c. 


T(y,0)  =  To(y) 


B.C. 


T(-»,t)  ■  T 


cs 


,2nrnn/’ 


ptu(p'Vm  -  q  ^  ) 

*  ^surface' 


kTyfO,*)  ~  1Burfaoe  r  + 


where 


r 


sur 
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There  are  several  differences  with  the  present  analysis,  but  the  functional 


form  of  the  unsteady  heat  transfer  boundary  condition  to  the  solid  is 
quite  similar.  This  system  was  solved  numerically  for  assumed  pressure 
variations  in  the  flame  which  simultate  ignition,  depressurization, 
and  oscillatory  motion.  The  results  show  the  possiblity  of  a  run-away 
burning  rate  under  certain  conditions  when  the  surface  reaction  is 
exothermic.  They  also  indicate  nonlinear  wave  forms  for  burning  rate 
as  a  function  of  time  for  sinusoidal  pressure  oscillations  in  the  flame. 
These  results  and  a  conclusion  that  the  extent  of  exothermic  reactions 
on  the  burning  surface  is  a  major  factor  in  determining  unstable  burning 
rates  are  supported  qualitatively  by  the  present  analysis  (see  Chapter 
IV). 

Also  discussed  in  Reference  (56)  is  the  important  point  con¬ 
cerning  the  effect  of  pressure  variations  in  the  flame  on  the  flame 
temperature.  It  is  often  argued  that  high  frequency  disturbances  occur 
adiabatically,  which  allows  the  instantaneous  flame  temperature  to  be 
computed  as  an  isentropic  function  of  pressure.  At  low  frequencies, 
the  outer  region  of  the  flame  is  usually  assumed  to  be  isothermal.  It 
is  quite  possible  that  neither  limiting  case  is  true  in  the  actual 
flame,  as  Summerfield  '  has  pointed  out.  Reference  (57)  speculates 

that  entropy  waves  may  appear  in  the  chamber  as  a  result  of  the  inter¬ 
action  between  pressure  disturbances  and  the  fluctuating  flame  temper¬ 
ature.  Krier,  et.  al.  were  the  first  to  present  estimates  of  the 

magnitude  of  this  entropy  variation  as  a  function  of  frequency.  Their 
results  (Figure  11,  Reference  56)  confirm  that  the  edge  of  the  flame 

zone  is  isothermal  under  conditions  of  a  low  frequency  pressure 
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oscillation.  Since  the  flame  temperature,  and  hence  the  enthalpy  of 
combustion,  hc,  is  assumed  to  be  constant  in  the  present  analysis,  the 
fundamental  chamber  frequency  was  computed  for  all  the  example  cases 
and  compared  to  Krier’s  predictions.  In  all  cases,  the  frequencies 
are  well  inside  the  range  where  the  isothermal  assumption  was  found 
to  be  valid. 

Krier's  combustion  model  is  extended  in  the  work  of  Merkle,  et. 

(59) 

al.  where  the  major  difference  is  the  treatment  of  the  gas -phase 

flame  region.  Reference  (59)  is  concerned  with  the  extinguishment  of 
burning  solid  propellants  and  hence  must  consider  the  burning  process 
at  pressures  much  lower  than  the  chamber  pressure.  The  improved  gas- 
ph'se  flame  analysis  is  based  on  Summerf ield ' s  Granular  Diffusion 

Flame  theory.  At  low  pressures  (and  low  temperatures)  the  flame  reaction 
is  assumed  to  be  kinetically  controlled  and  at  high  pressures,  it  is 
diffusion  controlled.  The  combination  of  these  mechanisms  allows  both 


processes  to  be  active  at  intermediate  pressure  levels,  but  the  limit¬ 
ing  cases  are  recovered  at  the  pressure  extremes.  The  analysis  also 
uses  an  Arrhenius -type  surface  reaction  law  in  place  of  the  previous 
power  law  temperature  dependence.  With  the  exception  of  some  assump¬ 
tions  regarding  the  gas-phase  reaction  rate,  this  analysis  and  the 
model  used  in  the  present  investigation  are  comparable  at  high  chamber 
pressures  where  diffusion  dominates.  The  results  presented  in  Reference 
(59)  describe  the  extinction  boundaries,  i.e.,  the  conditions  under 
which  the  combustion  process  is  quenched.  No  predictions  were  made 


for  oscillating  instability-type  pressure  disturbances. 


The  recent  study  of  Perctz  et.  al. 


presents  an  analysis  of 
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the  starting  transient  of  solid  propellant  rocket  motors  with  high 
internal  gas  velocities.  The  time-dependent  behavior  of  the  single 
phase  flow  in  the  combustion  chamber  is  found  with  a  numerical  solution 
of  the  one-dimensional  finite-difference  conservation  equations.  Pro- 

A 

visions  are  made  for  the  ignitor  gases  and  the  approximate  effect  of 
the  nozzle.  Below  a  given  value  of  ignition  temperature,  the  propellant 
grain  undergoes  a  heating  process  by  the  ignitor  gases,  described  by 


T  =  0 

yy 


I.C. 


T(y,0)  -  T0(y) 


B.C. 


T(-»,t)  -  Tcs 


kT  (0,t) 

y 


-  h  (T 
con  gas 


T 

sur 


) 


where  the  heat  transfer  coefficient,  hcQn,  is  an  empirical  expression 
accounting  for  many  effects.  The  problem  is  reduced  to  a  single 
ordinary  differential  equation  in  time  by  assuming  the  spatial  depen- 
dance  of  the  temperature  distribution  to  be  a  third-order  polynomial. 
The  resulting  ordinary  differential  equation  for  the  propellant  surface 
temperature  is  solved  simultaneously  with  the  chamber  flow.  When  the 
surface  temperature  exceeds  the  ignition  value,  the  burning  rate  is 


102 


assumed  to  be  predicted  by 


n  ..  -(^PSM 
r  =  ap  +  kh  e 
cor 

The  second  term  accounts  for  the  contribution  of  erosive  burning  if  the 
appropriate  constants  (k,hcon,3)  are  known.  The  final  results  for 
pressure  transients  in  the  engine  during  ignition  compare  well  with  the 
experimental  data  if  the  contribution  of  erosive  burning  is  included. 
Even  if  the  present  study  had  included  the  effect  of  erosive  burning 
in  the  above  manner,  comparisons  with  Reference  (60)  would  be  difficult 
due  to  the  steady  state  burning  expression  used  after  the  propellant  is 
ignited. 

Bie  related  field  of  armament  problems  also  requires  an  analysis 
of  the  solid  propellant  combustion  problem.  In  place  of  a  combustion 
chamber  lined  with  propellant  grain,  the  analysis  must  model  a  cylindri¬ 
cal  shell  partially  filled  with  a  granulated  propellant,  hence  referred 

(6lY 

to  as  a  porous  propellant.  In  Kuo'sv  ;  investigation  of  this  problem, 

several  complex  methods  of  predicting  the  time-dependent  heating  process 

of  the  porous  propellant  bed  are  discussed.  Since  the  analysis  is 

directed  toward  the  flame  spreading  phenomenon,  dynamic  burning  effects 

(55) 

are  ignored.  Thus,  similar  to  Most  ,  the  solution  to  the  time- 
dependent  heating  process  is  obtained  until  an  ignition  temperature  is 
achieved,  and  from  then  on,  the  burning  rate  is  predicted  with  the 
steady  state  expression,  r  =  apn.  The  time-dependent  problem  is  solved 
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with  an  integral  approach  which  assumes  the  temperature  distribution 
in  a  single  spherical  grain  is  given  by 


*(t)  ~ 

T(t,P.)  =  To  {e  P  “  Y(t)  T  } 

where 


=  propellant  grain  radius 


Substituting  this  expression  into  the  initial  value  -  boundary  value 
problem, 


-  !(BT>BR 


=  0 


I.C. 


T(R,0)  =  Tq(R) 


B.C. 


T(0,t)  =  Tcs 


TR(Bp,t)  = 


con 


[T 


gas 


-  T  ] 
sur 


results  in  a  single  ordinary  differential  equation  in  time  for  the 
unknown  coefficient,  Y(t), 
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—  =  function  (Y(t),  T  ,  R  ,  etc.) 
dt  '  '  /5  gas’  p’  ' 


The  solution  to  this  equation  is  obtained  simultaneously  with  the  gas 
dynamic  solution  of  the  flow  in  the  cylindrical  shell  which  provides 
T^s  •  The  present  study  was  unsuccessful  in  an  attempt  to  use  a 
similar  approach  because  of  the  complex  behavior  of  the  temperature 
distribution  in  the  solid  created  by  an  oscillatory  pressure  in  the 
flame. 

(62 ) 

Kitchens  '  '  has  studied  a  similar  problem  in  porous  propellants 

using  Kuo's  method  to  describe  the  heating  process,  and  hence  pre¬ 

dict  the  flame  spreading  process.  He  correctly  noted  that  for  a  spher¬ 
ical  propellant  grain  with  a  small  radius,  the  interior  boundary  con¬ 
dition  T(0,t)  =  T  ,  should  be  replaced  by  a  vanishing  radial  derivative, 
cs 

Chapter  IV  discusses  the  numerical  results  obtained  with  the 
combustion  model  derived  in  this  chapter,  as  well  as  the  chamber  and 
nozzle  flow  field  model  discussed  in  Chapter  II.  By  combining  these 
two  analyses,  predictions  of  the  transient  response  of  the  complete 
rocket  engine  with  a  burning  propellant  are  obtained. 


105 


CHAPTER  IV 


DISCUSSION  OF  RESULTS 

Numerical  predictions  of  the  response  of  a  rocket  engine  flow 
field  to  an  arbitrary  disturbance  are  obtained  with  a  combination  of 
two  computer  codes.  One  code  incorporates  the  combustion  chamber  and 
nozzle  flow  field  analysis  derived  in  Chapter  II.  The  other  code 
determines  the  time-dependent  surface  mass  flow  rate  at  a  given  burn¬ 
ing  station  along  the  solid  propellant  as  derived  in  Chapter  III.  The 
chamber  and  nozzle  flow  field  program  includes  the  option  to  accept  a 
previously-computed  flow  field,  superpose  a  continuous  or  shock  wave 
disturbance  on  this  flow  field,  and  then  restart  the  time  integration 
as  a  new  initial-value  problem.  In  the  examples  of  this  chapter,  the 
disturbances  are  introduced  on  the  steady  state  flow  field  for  the 
engine  under  consideration.  The  steady  state  flow  field  is  computed 
as  the  asymptotic  limit  of  the  time  integration,  starting  from  the 
assumption  that  the  propellant  is  ignited  and  burning  in  equilibrium 
at  ambient  pressure.  The  code  which  computes  the  time-dependent  burn¬ 
ing  rate  of  the  solid  propellant  can  be  separated  from  the  chamber 
and  nozzle  flow  field  program  and  applied  to  an  isolated  burning 
station.  The  combustion  response  curves  characteristic  of  each  pro¬ 
pellant  are  computed  with  this  code  by  imposing  a  particular  pressure 
oscillation  at  the  flame. 

The  full  transient  behavior  of  the  rocket  engine  flow  field  is 
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determined  with  a  combination  of  both  codes.  At  every  time  step,  the 
combustion  response  program  determines  the  propellant  mass  flow  rate 
at  each  of  the  burning  stations  along  the  grain.  The  boundary  conditions 
needed  for  determining  the  local  burning  rate  are  based  on  the  instan¬ 
taneous  chamber  flow  field  properties  directly  below,  i.e.,  at  the  same 
axial  location.  The  chamber  and  nozzle  flow  field  variables  are  then 
integrated  forward  in  time,  and  all  information  passing  through  the 
control  volume  perimeter  boundary  is  updated.  Thus,  the  solid  propellant 
combustion  process  is  solved  simultaneously  with  the  chamber  and  nozzle 
flow  field. 

One  of  the  complexities  of  the  flow  field  analysis  of  Chapter  II 
is  the  tracking  of  moving  shock  waves.  The  computer  routines  which 
perform  this  operation  were  checked  in  a  simple  test  case  describing 
the  flow  in  a  shock  tube.  In  this  computation,  the  effects  of  com¬ 
bustion,  solid  particles  and  area  variation  are  ignored.  A  constant 
area  chamber  with  closed  ends  is  divided  by  a  diaphram  separating 
regions  of  different  pressure.  At  t  =  0,  the  diaphram  is  removed  and 
the  time-dependent  calculation  tracks  a  shock  wave  propagating  into 
the  low  pressure  region.  The  flow  field  before  and  after  reflection 
from  the  closed  end  is  illustrated  in  Figure  14.  The  pressure  behind 
the  shock  wave  at  reflection  agrees  with  the  value  specified  in  a 
shock- tube  manual  '  This  computation  used  22  grid  points  and 

required  30  seconds  on  the  UNIVAC  1108.  *' 

Another  indication  of  proper  operation  of  the  computer  program 

is  its  ability  to  compute  a  steady  state  flow  field  as  the  asymptotic 

limit  of  a  time -dependent  computation.  For  the  arbitrarily  selected 
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engine  parameters  shown  as  Case  (l)  in  Appendix  B,  an  ignition  transient 
computation  was  started  with  the  assumption  that  the  pro,,  ellant  was 
burning  in  equilibrium  at  ambient  pressure.  The  propellant  mass  flow 
rate  from  the  surface  was  given  by 


m  =  0 .  Ip 
g 


0.3 


(4.1) 


instead  of  the  time  dependent  prediction  discussed  in  Chapter  III.  The 
pressure  at  the  head  end  of  the  motor  is  plotted  as  a  function  of  time 
in  Figure  15.  With  the  simple  burning  rate  aw,  apn,  there  is  no 
indication  of  oscillatory  motion  as  the  chamber  fills  and  the  pressure 
rises  monotonically  to  its  steady  state  value.  An  exception  was  found 
with  an  unrealistically  large  value  of  the  constant  "a".  This  produced 
a  very  large  initial  filling  rate  which  led  to  the  formation  of  a 
shock  wave;  the  shock  passed  through  the  throat  and  out  the  exit  plane, 
leaving  choked  flow  in  the  nozzle.  A  small  amplitude  wave  remained  in 
the  chamber  but  quickly  decayed.  The  results  in  Figure  15  confirm  that 
a  steady  flow  field  is  the  asymptotic  limit  of  the  time-dependent  inte¬ 
gration  in  the  computer  program.  Several  computer  experiments  indicated 
that  the  same  steady  state  flow  field  can  be  predicted  by  a  computation 
whose  initial  condition  is  a  reasonable  approximation  to  the  final  flow 
field.  This,  of  course,  would  be  expected.  Since  the  present  study  is 
directed  toward  the  transient  behavior  of  the  engine  after  it  has  achieved 
full  thrust,  computations  of  the  steady  state  flow  field  for  all  the 
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cases  to  follow  were  started  from  an  isentropic  flow  initial  condition. 
The  errors  in  this  approximate  initial  flow  field  were  quickly  eliminated 
in  the  time -dependent  computation.  All  flow  field  variables  of  the 
steady  state  solution  are  printed  on  cards  which  can  be  read  back  into 
the  program  at  a  later  time  to  restart  the  computation. 

The  steady  state  engine  flow  field  obtained  in  the  above  cal¬ 
culation  was  subjected  to  a  shock  wave  disturbance  to  study  the  engine 
response  characteristics.  Of  course,  the  new  calculation  has  exactly 
the  same  engine  parameters  and  combustion  model  as  the  previous  com¬ 
putation.  The  steady  state  flow  field  was  altered  with  the  addition 
of  a  shock  wave  disturbance  near  the  head  end  of  the  chamber,  simulating 
the  effect  of  an  explosive  charge  (see  t  =  0  flow  field  in  Figure  l6). 

The  pressures  behind  the  shock  wave  are  1.4  times  the  steady  state 
values  and  the  shock  begins  to  propagate  with  a  Mach  number  of  1.17. 
Figure  17  shows  the  pressure-time  history  at  the  head  end  of  the  chamber. 
The  disturbance  is  still  a  shock  wave  on  the  first  reflection  but  all 
other  reflections  indicate  a  continuous  disturbance.  Spatial  pressure 
distributions  in  the  motor  at  the  four  time  points  indicated  in  Figure 
17  are  illustrated  in  Figure  16.  The  mass  flow  control  of  a  choked 
throat  is  shown  by  the  pressure  profile  at  point  (l)  .  The  shock  wave 
has  passed  through  the  throat  leaving  behind  a  large  surge  of  pressure 
which  is  effectively  blocked  by  the  choked  condition.  At  this  instant, 
the  Mach  number  in  the  geometric  throat  is  approximately  1.2,  as  shown 
in  Figure  18.  This  would  imply  that  an  approximate  nozzle  entrance- 
plane  boundary  condition  which  assumes  a  unity  throat  Mach  number  intro¬ 
duces  unknown  error  into  the  reflection  process.  As  the  pressure  surge 
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(in  ©)  propagates  back  upstream,  it  reforms  into  a  shock  wave  which 
reflects  from  the  head  end  of  the  chamber.  The  strength  of  the  reflected 
shock  is  greatly  reduced,  as  shown  in  ©  .  After  the  second  nozzle 
reflection  process,  the  remaining  disturbance  is  too  weak  to  reform 
into  a  shock  wave  as  shown  in  ©  .  The  disturbance  remains  a  continuous 
wave  from  here  on.  The  distribution  at  point  ©  also  indicates  the 
appearance  of  a  higher  mode  which  is  following  the  main  disturbance. 

This  is  also  confirmed  by  the  time-history  in  Figure  17  with  the  smaller 
pressure  wavelet  (at  t  =  1.4)  which  is  following  the  main  disturbance. 

At  time  point  ©  ,  the  pressure  distribution  shows  that  only  small 
amplitude  disturbances  remain.  The  rapid  decay  of  the  original  shock 
wave  disturbance  is  due  to  two  effects.  The  convective  losses  through 
the  nozzle  are  very  large  since  the  area  ratio  is  4.0,  which  is  unreal¬ 
istically  low  for  most  solid  proA  pliant  rocket  engines .  In  addition, 

the  combustion  response  of  the  propellant  was  predicted  with  the  stc  iy 

0  3 

state  expression,  m  =  O.lp  "  .  The  ability  of  this  model  to  drive  or 

S 

sustain  a  disturbance  traveling  in  the  combustion  chamber  is  minimal. 

This  point  has  been  conjectured  before,  and  it  is  borne  out  by  the 
present  computation. 

The  effect  of  solid  particles  in  the  gas  flow  on  the  engine 
response  was  investigated  in  a  similar  fashion.  The  steady  state  com¬ 
bustion  expression,  apn,  was  used  with  a  pressure  exponent  of  0.4. 

The  area  ratio  of  the  nozzle  was  i' creased  to  11.6  to  reduce  the  large 
convective  losses  of  the  previous  case.  The  behavior  of  two  nearly 
identical  engines  (see  Case  2,  Appendix  B)  was  investigated,  one  with¬ 
out  particles  and  the  other  with  10|i  average  diameter  particles  at  a 
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ten  percent  weight  loading.  A  continuous  disturbance  (see  Figure  19), 
whose  maximum  amplitude  was  forty  percent  of  the  chamber  pressure  was 
superposed  on  the  steady  state  pressure  distribution  in  both  cases. 

The  head  end  pres  sure -time  history  of  the  altered  engine  flow  field 
with  no  particles  is  shown  in  Figure  20.  The  decay  rate  of  the  dis¬ 
turbance  is  less  than  in  Case  (l)  due  to  the  increased  area  ratio, 
but  the  combustion  response  is  again  too  weak  to  sustain  the  distur¬ 
bance.  The  head  end  pressure-time  history  of  the  similar  engine  with 
10p  diameter  particles  is  shown  in  Figure  21.  The  dashed  curve  in 
this  figure  is  the  same  result  shown  in  Figure  20.  Comparison  of  the 
two  curves  reveals  several  anticipated  trends  resulting  from  the 
addition  of  solid  particles  to  the  flow.  In  this  example,  the  small 
amount  of  higher  harmonic  content  when  no  particles  are  present  is 
smoothed  out  when  the  10p,  particles  are  added  to  the  flow .  Thus, 
disturbances  at  higher  frequencies  than  the  first  fundamental  mode 
may  be  strongly  damped  by  the  drag  effect  of  the  particles,  if  the 
average  particle  diameter  is  tuned  to  that  frequency.  Although  this 
smoothing  effect  broadens  the  base  of  the  pressure  wave-form  at  the 
fundamental  frequency,  it  does  not  noticeably  affect  the  decay  rate 
of  the  pressure  maximums.  Another  prominant  feature  is  the  decrease 
in  the  effective  wave  speed  in  the  gas-particle  mixture  over  a  cycle, 
as  demonstrated  by  the  increased  time  between  reflections  at  the 
head  end. 

The  computer  results  in  the  above  example  confirm  the  expected 
effect  of  solid  particles.  They  also  serve  to  emphasize  that  modeling 
these  particles  as  inert  spheres  which  create  drag  and  dissipate  energy 
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will  not  provide  an  explanation  of  why  longitudinal  instabilities  are 
sometimes  aggrevated  with  the  use  of  metalized  solid  propellants. 
Although  spherical  inert  particles  can  be  expected  to  attenuate  the 
high  frequency  unsteady  motion  and  contribute  to  the  broadening  of  lower 
frequency  large  amplitude  motion,  their  overall  importance  in  the 
present  problem  is  secondary.  For  this  reason,  further  results  des¬ 
cribed  in  this  report  will  neglect  the  presence  of  solid  particles 
in  the  flow.  Consequently,  any  conclusions  drawn  from  results  snowing 
a  decaying  disturbance  without  particles  present  remain  valid;  however, 
conclusions  about  the  growth  of  disturbances  must  be  tempered  with  the 
possible  contribution  of  the  particle  flow. 

The  behavior  of  the  combustion  response  was  investigated  for 
an  isolated  burning  station  using  the  time-dependent  model  derived  in 
Chapter  III.  In  each  of  the  following  examples,  a  given  pressure-time 
history  was  imposed  on  the  gas-phase  flame.  Hence,  it  is  assumed  that 
the  mass  flow  response  from  the  surface  has  no  effect  on  the  pressure 
disturbance;  this  is  not  the  case  when  the  burning  stations  are  coupled 
to  the  chamber  and  nozzle  flow  field.  For  the  set  of  propellant  ,and 
flame  parameters  listed  as  Case  (3)  in  Appendix  B,  the  time -dependent 
surface  mass  flow  rate  was  computed  for  two  sinusoidal  pressure  oscil¬ 
lations,  each  with  an  amplitude  of  ten  percent  of  the  pressure  in  the 
flame.  This  example  demonstrates  the  difference  in  combustion  response 
to  a  xow  frequency  (w  =  2tt)  and  a  high  frequency  (w  =  ^8tt)  disturbance. 
Using  a  time-scale  non-dimensionalized  by  the  period  of  oscillation, 
the  results  in  Figure  22  show  that  the  mass  flow  rate  "leads"  the  low 
frequency  pressure  oscillation  and  "lags"  the  high  frequency  one.  This 
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is  in  agreement  with  the  results  from  linearized  combustion  response 
theory  To  demonstrate  some  of  the  nonlinear  character  of  the 

time-dependent  model,  the  two  cases  were  rerun  with  a  ten  percent 
amplitude  square  wave  pressure  disturbance  of  the  same  period.  The 
response  to  the  square  wave  and  the  sine  wave  are  shown  together  in 
Figure  23  for  the  low  frequency  case  and  in  Figure  24  for  the  high 
frequency  case.  The  responses  are  distinctly  nonlinear.  These  results 
support  a  conclusion  from  Reference  (56)  that  the  nonlinear  effect  on 
peak  to  peak  mass  flow  rate  variations  is  diminished  at  high  frequency. 
The  square  wave  pressure  disturbance  is  admittedly  artificial,  but  the 
results  demonstrate  an  important  point.  When  the  pressure  disturbance 
is  no  longer  a  small  amplitude  sinusoidal  oscillation,  the  nonlinearity 
of  the  system  (Equation  3-20)  is  very  important  in  determining  the 
transient  combustion  response.  Thus,  predicting  the  response  to 
arbitrary  pressure  disturbances  with  the  results  of  a  linearized  theory 
could  lead  to  large  errors. 

One  of  the  important  results  from  linearized  combustion  response 

theory  is  a  plot  of  the  real  part  of  the  response  function,  R  E  ? 

Ap/p 

as  a  function  of  the  nondimensionalized  frequency,  H.  A  curve  of  this 
type  can  be  constructed  for  each  propellant  system.  The  value  at  a 
given  frequency  represents  the  magnitude  of  the  response  to  a  small 
amplitude  pressure  oscillation,  but  it  is  not  an  indication  of  whether 
a  rocket  engine  containing  this  propellant  will  be  stable  or  unstable. 
It  is  a  useful  guide  to  estimating  the  effect  of  a  change  in  mean 
burning  rate,  chamber  temperature,  etc.,  on  the  behavior  of  the  pro¬ 
pellant  in  the  engine.  A  similar  combustion  response  curve  can  be 
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predicted  with  the  nonlinear  timer dependent  model  by  repeating  the 
computations,  shown  in  Figure  22  for  small  amplitude  pressure  oscillations, 
over  a  wide  frequency  range.  A  value  of  ft  is  computed  from  each  time- 
dependent  profile.  The  symbols  in  Figure  25  are  the  results  of  such 
a  calculation.  Also  plotted  on  the  graph  is  the  two-parameter  combustion 
response  function  predicted  by  linear  theory  for  the  combustion  model 
used  in  the  time -dependent  solution.  The  values  of  A  and  B  were  computed 
from  Equations  (4.l4)  and  (4.32)  of  Reference  (34).  The  agreement  is 
good  except  at  the  higher  frequencies.  However,  for  other  propellant 
systems  to  be  discussed,  a  similar  difference  between  the  two  computations 
also  appears  at  low  frequencies.  It  is  believed  that  these  discrepancies 
are  due  to  the  influence  of  the  time-dependent  regression  rate  multi¬ 
plying  the  temperature  gradient  term  in  Equation  (3 .20)  and  the  behavior 
of  the  nonlinear  boundary  condition,  Equation  (3.20)(c). 

It  should  be  emphasized  that  not  all  propellant  systems  respond 
to  a  ten  percent .  amplitude  sinusoidal  pressure  oscillation  with  a 
sinusoidal  mass  flow  rate.  The  propellant  parameters  used  in  the  cal¬ 
culation  of  Figure  25  lead  to  what  might  be  called  a  "mild"  response, 
since  the  value  of  the  maximum  response  (where  the  mass  flow  rate  is 
in-phase  with  the  pressure  oscillation)  is  about  2.2.  Other  more 
responsive  propellant  systems  may  exhibit  an  entirely  different  com¬ 
bustion  response,  as  will  be  evident' shortly. 

Predictions  of  engine  stability  require  the  simultaneous  solution 
of  the  unsteady  combustion  process  and  the  unsteady  flow  in  the  chamber 
and  nozzle.  A  large  number  of  propellant  and  engine  parameters  must 
be  specified  in  any  one  case.  With  the  limited  time  available,  the 
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present  analysis  chose  to  study  the  influence  of  the  pyrolysis  reaction 
at  the  propellant  surface  upon  engine  stability.  There  is  considerable 
speculation  about  the  nature  and  extent  of  this  reaction.  Beckstead, 

Derr  and  Price  have  investigated  the  steady  state  burning  rate  of 

monopropellants  and  composite  propellants  with  AP  as  an  oxidizer.  Their 
results  suggest  that  the  net  heat  release  of  this  reaction  is  exothermic, 
with  a  value  of  approximately  120  cal/gm  for  an  AP  monopropellant.  Since 
the  heat  of  sublimation  for  AP  is  about  585  cal/gm  (endothermic),  a  very 
exothermic  reaction  is  occurring  at  the.  surface.  In  an  attempt  to 
isolate  the  effect  of  this  reaction,  the  following  combustion  parameters 
were  held  constant  in  the  examples  which  follow: 

Solid  propellant: 

p  a  density  of  propellant  =  1.95  gm/cnr  =  121.8  lbm/ft3 
s 

Tcs  -  cold  solid  temperature  =  300° K  =  54o°R 

E  =  activation  energy  of  surface  reaction  =  20  kcal/mole 
s 

-X* 

cg  =  heat  capacity  of  propellant  =  .275  cal/gm°K  =  .275  BTU/lbm-°R 

*  -3/o 

k  =  thermal  conductivity  of  propellant  =  1.20  x  10  cal/cm-sec-  K 

s 

=  8.06  x  10"5  BTU/ft-sec-°R 

<• 

Gas  phase: 

V  =  isentropic  index  =  1.2 

0^  =  specific  heat  =  0.300  cal/gm-°K  =  0.300  BTU/lbm-°R 

T^,  =  flame  temperature  =  2500°K  =  4500°R 

W*  =  average  molecular  weight  =  22  gm/gm-mole  =  22  lb^/lb  -mole 

E^,  s=  activation  energy  of  flame  reaction  =  30  kcal/mole 

*  -4/o 

^  =  thermal  conductivity  =  2.0  x  10  cal/cm-sec-  K 
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The  remaining  parameters  are  determined  by  balancing  the  equations  at 

steady  state  with  typical  experimental  values.  For  example,  a  regression 

rate  of  1  cm/sec  at  a  surface  temperature  of  880°K  determines  the  con- 

stant  Bf  in  Equation  (3»3).  Using  this  surface  temperature  and  the 

* 

assumed  value  for  Q  the  steady  state  energy  balance  across  the  pro- 

S  “ 

* 

pellant  and  the  flame  specifies  the  heat  release  term,  Q^.,  in  the  gas 

phase.  Then  with  the  value  of  the  expected  chamber  pressure  in  the 

engine  at  steady  state,  the  heat  transfer  boundary  condition  at  the 

propellant  surface  determines  the  constant  in  the  gas  phase  reaction 

rate,  i.e.,  uj  (see  Equation  3.1l)»  A  sample  calculation  is  shown 
con 

in  Appendix  C. 

A  starting  point  in  the  investigation  of  nonlinear  stability  of 

the  complete  rocket  engine  was  a  propellant  system  with  only  sublimation 

as  a  surface  reaction.  Using  the  above-mentioned  propellant  parameters 

* 

which  are  typical  of  AP  propellants,  and  the  assumption  that  Q  =  +  5o5 

s 

cal/gm  (+  1,053  BTU/lb^)  endothermic,  the  steady  state  balance  procedure 
was  used  to  determine  the  remaining  parameters  [see  Case  (4),  Appendix 

B,  also  Appendix  C],  These  values  were  used  to  compute  A  and  B  for 

* 

the  two-parameter  combustion  response  function  which  is  plotted  in 
Figure  26.  The  curve  is  quite  flat,  indicating  no  preferential  range 
of  frequencies.  In  this  case,  a  small  amplitude  oscillation  at  the 
fundamental  frequency  of  the  chamber,  corresponding  to  0  =  10  in  Figure 
26,  produces  a  value  of  mass  flow  response  less  than  0.9*  The  steady 
state  flow  field  for  this  engine  was  altered  with  a  shock  wave  distur¬ 
bance  in  the  same  manner  as  the  t  =  0  plot  of  Figure  16.  The  pressure 

increase  behind  the  shock  is  forty  percent  of  the  steady  state  pressure 
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at  the  same  location.  The  pressure-time  history  at  the  head  end  of  the 
chamber  is  shown  in  Figure  27.  The  results  show  that  the  initial  shock 
wave  has  become  a  continuous  disturbance  on  the  first  reflection,  and 
continues  to  decay  thereafter.  The  effect  of  the  flat  combustion  res¬ 
ponse  curve  is  evident  as  the  decaying  fundamental  mode  transfers 
energy  to  the  higher  harmonics.  Since  all  modes  have  nearly  the  same 
value  of  combustion  response,  the  higher  harmonics  remain  pi'ominent  as 
the  wave  form  decays .  It  is  concluded  that  this  engine  is  stable  when 
perturbed  by  a  large  amplitude  disturbance  and,  therefore  it  is  unlikely 
that  this  system  would  be  capable  of  spontaneous  instability  in  the 
longitudinal  mode. 

The  effect  of  an  exothermic  surface  reaction  on  the  stability 
of  a  rcckct  engine  is  studied  in  the  following  example.  The  combustion 
response  curve  in  Figure  25  is  the  result  of  assuming  the  net  heat 
release  in  the  pyrolysis  reaction  to  be  100  cal/gm,  exothermic.  Two 
engines  were  designed  so  that  their  fundamental  frequencies  operate 
at  different  points  on  the  same  curve.  Referring  to  Figure  25,  Case 
(5)  operates  at  0  =  10  with  a  response  of  approximately  2.  Since  this 
point  is  past  the  peak  response,  the  propellant  mass  flow  rate  would 
be  expected  to  lag  a  small  amplitude  pressure  oscillation  at  the 
fundamental  frequency  (and  for  higher  harmonics).  Case  (6)  operates 
at  fi  =  1.5  with  a  response  of  1.4.  At  this  point,  the  propellant  mass 
flow  rate  would  be  expected  to  lead  a  small  amplitude  pressure  oscillation 
at  the  fundamental  frequency.  It  would  also  lead  at  the  next  two  har¬ 
monics,  each  of  which  has  a  larger  value  of  combustion  response  than  - 
the  fundamental  mode.  Both  engines  were  subjected  to  an  initial  shock 
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wave  disturbance  which  created  a  forty  percent  increase  in  pressure 
(similar  to  t  =  0  plot  of  Figure  l6).  The  head  end  pressure-time 
history,  along  with  the  instantaneous  mass  flow  response  for  the  burning 
station  located  at  the  head  end,  are  shown  together  in  Figure  28  for 
Case  (5)  and  in  Figure  29  for  Case  (6).  Comparison  of  the  two  figures 
indicates  a  noticeable  difference  in  the  mass  flow  response.  For  a 
large  amplitude  disturbance  (particularly  a  shock  wave),  the  concept 
of  leading  or  lagging  the  pressure  oscillation  is  augmented  by  the 
difference  in  the  transient  response's.  Case  5  (lag)  responds  with 
wider  and  more  gentle  mass  flow  rate  peaks,  and  with  very  little  higher 
harmonic  content  as  the  pressure  wave  form  decays.  The  minimum  amount 
of  higher  mode  content  in  the  pressure  wave  form  is  the  result  of  the 
reduced  mass  flow  response  at  higher  frequencies  with  an  even  greater 
amount  of  lag.  Case  (6)  (lead)  shows  a  sharp-peaked  mass  flow  rate 
which  responds  right  behind  the  shock  wave,  further  increasing  the 
pressure  at  reflection.  This  acts  as  a  driver  for  the  shock  wave 
which  is  sustained  for  two  cycles.  Once  the  initial  shock  becomes 
a  large  amplitude  continuous  wave,  the  mass  flow  response  does  not 
lead  it,  but  remains  nearly  in-phase  with  it.  In  addition  to  the 
greater  magnitude  of  response  at  higher  frequencies,  this  creates  a 
more  prominent  higher  harmonic  content  in  the  pressure  waveform.  The 
fact  that  the  peak  value  of  the  disturbance  after  six  cycles  is  slightly 
greater  in  Case  (6)  [ft  =  1.4]  than  Case  (5)  [ft  =  2.0]  demonstrates  the 
importance  of  the  phase  condition.  In  general,  a  mass  flow  rate  which 
responds  after  the  pressure  disturbance  has  passed  is  not  an  effective 


driver  for  combustion  instabilities. 
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The  rocket  engines  analyzed  so  far  have  been  unable  to  sustain 
large  amplitude  disturbances  and  would  be  considered  quite  stable.  For 
the  engine  parameters  used  in  these  examples,  the  combustion  response 
is  insufficient  to  match  the  losses  through  the  nozzle.  If  the  exo¬ 
thermic  heat  release  in  the  surface  reaction  is  adjusted  to  115  cal/gm, 
the  response  function  increases  to  that  shown  in  Figure  30.  The 
following  example  [Case  (7),  Appendix  B]  employs  this  propellant  system 
in  the  rocket  engine  of  Case  (5).  Again,  the  operating  point  for  the 
fundamental  mode  is  n  =  10,  but  the  combustion  response  magnitude  has 
increased  to  approximately  3*  as  computed  by  the  nonlinear  time -dependent 
combustion  model.  The  steady  state  flow  field  for  this  motor  was  altered 
with  a  continuous  disturbance  with  a  maximum  amplitude  of  forty  percent 
of  the  chamber  pressure  (same  as  Figure  19).  The  head  end  pressure¬ 
time  history  of  the  restarted  computation  is  illustrated  in  Figure  31. 

At  the  end  of  the  first  cycle  the  amplitude  of  the  disturbance  has 
increase!,  but  it  then  proceeds  to  decay.  Thus  even  this  response 
function  is  not  sufficient  to  sustain  the  disturbance  when  the  mass 
flow  response  lags  the  pressure  disturbances  at  all  chamber  frequencies. 

Increasing  the  net  heat  release  of  the  surface  pyrolysis  reaction 
from  115  cal/gra  to  122  cal/gm  (exothermic)  produces  a  large  increase 
in  the  combustion  response  curve,  as  shown  in  Figure  32.  This  combustion 
model  was  used  in  the  same  rocket  engine  as  discussed  in  the  previous 
paragraph.  The  nonlinear  time -dependent  solution  for  the  combustion 
response  indicates  a  response  value  of  approximately  3.6  at  the  funda¬ 
mental  frequency  (f)  =  10 )  of  the  combustion  chamber  with  the  new  model  - 
[see  Case  (8),  Appendix  B].  The  steady  state  flow  field  was  altered 
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with  a  continuous  disturbance  with  a  maximum  amplitude  of  twenty  per¬ 
cent  of  the  chamber  pressure  (see  Figure  33).  The  head  end  pressure¬ 
time  history  of  the  restarted  computation  is  shown  in  Figure  3^.  The 
disturbance  increases  at  first,  and  then  decays  slowly.  The  out-of¬ 
phase  mass  flow  response  is  still  evident  in  the  behavior  of  the  higher 
modes.  To  investigate  the  possibility  of  triggering  this  engine  with 
a  larger  amplitude  disturbance,  the  computation  was  rerun  with  a  forty 
percent  amplitude  continuous  disturbance  (Figure  33)  of  the  same  shape. 
The  head  end  pressure-time  history  along  With  the  instantaneous  mass 
flow  rate  for  the  burning  station  at  the  same  location  are  shown  in 
Figure  35.  This  plot  indicates  the  possible  formation  of  a  limit  cycle 
oscillation  at  nearly  twice  the  initial  disturbance  amplitude.  However 
a  dramatic  change  is  evident  in  the  character  of  the  mass  flow  response 
This  behavior  is  a  sharp  but  finite-amplitude  "spike".  As  would  be 
expected,  when  the  burning  rate  spikes,  the  flow  field  pressure  also 
rises  sharply  at  the  same  location. 

The  phenomenon  of  a  burning  rate  spike  was  not  anticipated  when 
the  computational  method  was  constructed.  Although  provisions  were 
made  for  traveling  discontinuities,  the  region  between  any  two  moving 
boundaries  was  assumed  to  be  smooth  and  continuous.  Hence,  a  Taylor 
series  truncated  after  the  second  order  term  is  the  basis  for  the 
forward  marching  time  integration  cf  the  conservation  equations.  As 
an  illustration,  consider  the  homogeneous  equations  for  gas  flow  only: 


f ,  +  Af  =  0 
t  x 
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where  f  is  a  vector  with  components  p,  u,  and  S,  and  A  is  the  appropriate 
coefficient  matrix.  If  this  system  is  expanded  in  a  second-order  Taylor 
series,  the  truncation  "error"  is  given  by^^ 

\  Kxx  ^  *  fttt  it3] +  •  •  • 

If  a  mass  flow  spike  creates  a  similar  disturbance  in  f,  then  f  and 
F  ’  ttt 

f  can  become  significant.  Tlius,  integrating  the  flow  field  equations 

XXX 

with  the  present  scheme  introduces  unknown  truncation  errors  whenever 
the  mass  flow  spikes .  There  is  no  obvious  way  to  properly  account  for 
this  phenomenon  within  the  present  computational  method  without  recourse 
to  unreasonably  small  values  of  At  and  Ax.  Thus,  flow  field  predictions 
which  include  mass  flow  rate  spikes  are  assumed  to  represent  the  proper 
trend  in  the  rocket  engine,  but  the  numerical  values  have  unknown  accuracy. 
Hence,  no  attempt  was  made  to  compute  the  final  form  of  a  limit  cycle 
oscillation  in  the  motor  under  these  conditions. 

The  behavior  of  the  nonlinear  combustion  model  can  be  studied 
in  more  detail  without  the  combustion  chamber  analysis  and  its  asso¬ 
ciated  truncation  error.  The  following  results  were  generated  by  the 
computer  code  which  solves  the  time-dependent  burning  problem  at  a 
single  burning  station  for  any  assumed  pressure  variation  in  the  flame 
zone.  Eacli  calculation  used  1000  time  "steps"  per  cycle  of  pressure 
oscillation  and  the  convergence  criterion  for  satisfying  the  boundary 
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condition  at  the  interface  was  seven  significant  digits.  The  response 

of  the  combustion  model  used  in  Case  (8)  [i.e.,  Figure  32,  Q  =  -  122 

s 

cal/gm]  was  determined  for  several  sine  wave  pressure  oscillations  at 
the  fundamental  frequency  of  the  chamber  (fl  =  10).  Figure  36  shows 
the  mass  flow  response  to  a  sinusoidal  pressure  oscillation  with  an 
amplitude  of  five  percent  of  the  mean  chamber  pressure  (note:  the 
time  scale  is  nondimensionalized  by  the  period  of  the  pressure  oscillation, 
i.e.,  t  =  5  is  the  end  of  the  fifth  cycle).  After  a  transient  adjustment 
during  the  first  cycle,  the  mass  flow  rate  establishes  a  sinusoidal -type 
wave  form.  Closer  inspection  reveals  that  the  maximum  values  of  this 
wave  form  alternate  between  a  high  and  low  value  for  successive  cycles. 

The  minimum  values  are  the  same.  The  alternating  high-low  trend  is 
very  pronounced  when  the  amplitude  of  the  pressure  oscillation  is 
increased  to  ten  percent  as  shown  in  Figure  37.  This  result  shows 
a  transient  behavior  which  leads  to  a  spike  on  the  fifth  cycle.  Repeat¬ 
ing  the  same  calculation  with  a  fifteen  percent  amplitude  pressure 
oscillation  produces  an  alternating  pattern  of  spikes  on  every  other 
cycle  as  shown  in  Figure  38.  Finally,  increasing  the  amplitude  of  the 
pressure  oscillation  to  twenty  percent  results  in  a  similar  wave  form 
which  repeats  every  two  cycles  (see  Figure  39)- 

These  predictions  of  combustion  response  are  based  on  sinusoidal 
pressure  oscillations  at  a  single  frequency,  corresponding  to  n  =  10  on 
the  response  curve  shown  in  Figure  32.  Since  the  surface  mass  flow 
rate  lags  a  small  amplitude  pressure  oscillation  at  this  frequency, 
it  was  decided  to  investigate  the  lead  portion  of  the  curve  in  the 
same  way.  A  lead  case  was  chosen  at  ft  =  3  because  the  combustion 
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response  magnitude  is  almost  identical  to  that  at  n  -•  10.  The  pro¬ 
pellant  mass  flow  rate  due  to  a  five  percent  amplitude  sinusoidal 
pressure  oscillation  at  0  =  3  is  shown  in  Figure  40  (note:  the  time 
scale  is  nondimens ionali zed  by  the  period  of  the  pressure  oscillation). 
The  response  to  a  six  percent  pressure  oscillation  is  also  shown  on 
the  same  graph.  Comparison  of  the  two  curves  illustrates  the  evolution 
of  the  spike  and  the  small  amplitude  peak  seen  in  the  previous  results. 
However,  the  phase  characteristic  plays  an  important  role.  For  pressure 
oscillations  at  0  -  10,  the  spikes  formed  after  the  pressure  maximum; 
at  0  =  3>  they  appear  to  be  forming  before  the  pressure  maximum.  This 
is  confirmed  when  the  pressure  amplitude  is  increased  to  ten  percent 
as  shown  in  Figure  4l.  Note  that  this  wave-form  is  similar  to  that  in 
Figure  39  which  required  a  twenty  percent  amplitude  pressure  oscillation. 
However,  the  oscillation  at  D  =  3  leads  to  a  spike  on  every  cycle, 
where  at  0  =  10  the  spike  occurred  on  every  other  cycle. 

It  is  also  of  interest  to  determine  the  nonlinear  behavior  of 

the  combustion  model  at  the  frequency  corresponding  to  the  maximum 

response  value.  This  was  carried  out  for  a  combustion  model  with  a 

net  heat  release  at  the  surface  of  120  cal/gra  (exothermic).  The 

corresponding  combustion  response  curve  js  shown  in  Figure  42,  which 

indicates  that  the  mass  flow  rate  should  be  in-phase  with  a  small 

amplitude  pressure  oscillation  at  a  frequency  of  Q  =  5.5.  The  mass 

flow  response  due  to  a  sinusoidal  pressure  oscillation  at  Q  =  5.r> 

with  an  amplitude  of  one  percent  of  the  mean  pressure  in  the  flame  is 

shown  in  Figure  43a.  After  a  short  transient  period,  the  wave  form 

shows  no  evidence  of  the  nonlinear  influence.  Repeating  this 
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calculation  with  a  five  percent  amplitude  pressure  oscillation  produces 
the  distorted  wave-form  shown  in  Figure  43b.  This  closely  resembles 
the  result  shown  in  Reference  (5 6)  for  a  ten  percent  pressure  oscillation 
in  a  combustion  model  with  H  =  .80  (H  =  -  Q  /c  (T  -  T  )  =  0.75  in  the 
present  example).  If  the  pressure  amplitude  is  increased  to  7.5  percent, 
the  mass  flow  response  behaves  in  the  transient  manner  shown  in  Figure 
43c.  A  spike  occurs  just  after  the  pressure  maximum  on  the  second  cycle, 
but  does  not  reappear  until  the  fifth  cycle.  This  behavior  seems  to  indi¬ 
cate  the  existence  of  a  threshold  amplitude  criterion,  possibly  including 
a  time  derivative,  which  must  be  exceeded  before  a  spike  can  occur  on 
every  cycle.  Whatever  the  criterion,  it  is  exceeded  when  the  pressure 
oscillation  amplitude  is  increased  to  ten  percent,  as  shown  in  Figure 
43d.  After  the  transient  adjustment  of  the  first  cycle,  the  mass  flow 
response  is  characterized  by  a  repeating  pattern  of  spikes.  In  the 
response  to  a  7*5  percent  amplitude  oscillation,  the  two  mass  flow  rate 
spikes  occurred  after  the  pressure  maximums.  In  the  response  to  a  ten 
percent  amplitude  oscillation,  the  spikes  are  exactly  in-phase  with 
the  pressure  maximums.  This  type  of  combustion  response  would  exert 
a  substantial  driving  force  on  a  disturbance  in  the  combustion  chamber. 

The  combustion  response  spike  can  be  examined  in  more  detail  with 
a  semi-log  plot  on  an  expanded  time  scale.  The  mass  flow  rate  between 
t  =  3.20  to  3-30  shown  in  Figure  43d  is  replotted  in  Figure  44.  This 
time  interval  corresponds  to  0.543  mil-sec  in  real  time.  Most  of  the 

4*  Vi  +  Vi 

"spike"  occurs  in  l/lCr  of  this  time,  of  l/KHj  of  the  pressure 
oscillation.  Since  the  spike  takes  place  over  a  short  time  interval, 
it  would  be  expected  that  the  related  temperature  fluctuations  in  the 
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solid  propellant  will  be  restricted  to  a  region  close  to  the  surface. 

Three  temperature  profiles  in  the  unburned  solid  are  shown  in  Figure 
45  for  the  "spike"  illustrated  in  Figure  44.  The  temperature  profile 
©  at  t  =  3-022  is  before  the  peak,  profile  (?)  at  t  =  3.245  is  at  the 
peak,  and  profile  Q)  at  t  =  3.283  occurs  after  the  peak.  The  temper¬ 
ature  profile  at  t  =  4.022,  i.e.,  one  cycle  later,  is  identical  to  that 
in  Q)  .  During  the  spike,  the  movement  of  the  time-dependent  temperature 
distribution  in  the  solid  is  seen  to  be  confined  within  three  non- 

dimensional  unit-depths  from  the  surface.  For  this  propellant  system, 

_2 

three  units  represents  1.3  x  10  cm.  Actually,  the  large  amplitude 

temperatures  propagate  inward  less  than  l/lO  of  this  distance.  Thus, 

temperature  fluctuations  very  close  to  the  propellant  surface  play  a 

dominant  role  in  the  unsteady  mass  flow  response. 

The  temperature  distribution  at  the  peak  of  the  spike  shows  that 

a  severe  gradient  at  the  surface  is  sustained  for  a  short  time.  If 

the  thermal  conductivity  of  the  unburned  propellant  were  increased,  this 

should  produce  a  noticeable  effect  on  the  surface  temperature  fluctuations . 

The  propellant  system  (Figure  42)  which  exhibits  the  repeating  spikes 

shown  in  Figure  43d  was  altered  by  doubling  the  thermal  conductivity  of 

the  solid  propellant  (ko  =  2.4  x  10  cal/sec-  K-cm).  All  other 

°new 

physical  properties  remained  the  same.  This  new  system  responds  to 
the  same  ten  percent  amplitude  pressure  oscillation  in  the  flame  with 
the  v  as s  flow  rate  shown  in  Figure  46.  The  results  show  no  trace  of 
any  spike.  This  change  can  also  be  reasoned  from  the  combustion  res¬ 
ponse  curve  shown  in  Figure  42.  Doubling  the  thermal  conductivity  does 
not  change  the  response  curve,  but  it  does  shift  the  frequency  of  the 
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operating  point.  Thus,  the  results  in  Figure  46  are  for  a  pressure 

disturbance  at  0  =  11 »  which  is  well  into  the  lag  portion  of  the  curve. 

The  fact  that  Figure  46  resembles  Figure  36  is  to  be  expected.  Another 

important  parameter  governing  the  mass  flow  response  is  the  activation 

energy  of  the  surface  pyrolysis  reaction.  Decreasing  the  value  of  the 

activation  energy  will  substantially  reduce  the  sensitivity  of  the 

regression  rate  to  temperature  variations.  The  propellant  system  with 

the  response  curve  shown  in  Figure  42  was  again  altered  by  assuming 

the  activation  energy  of  the  surface  reaction,  E  ,  to  be  10  kcal/mole. 

This  new  system  was  subjected  to  the  same  ten  percent  amplitude  pressure 

oscillation  which  led  to  the  repeating  spikes  of  Figure  43d.  The  mass 

flow  response  is  shown  in  Figure  47  to  be  a  nearly  perfect  sine  wave. 

This  change  might  also  be  anticipated  from  the  movement  of  the  combustion 

response  curve.  Plotted  on  Figure  42  is  the  linear  combustion  response 

curve  for  the  above  propellant  system  with  Eg  =  10  kcal/mole.  The 

operating  frequency  is  still  0  =  5-5>  but  the  magnitude  of  the  response 

is  only  slightly  greater  than  unity.  This  accounts  for  the  sinusoidal 

mass  flow  response  shown  in  Figure  47. 

In  the  study  up  to  this  point,  all  computations  of  the  transient 

behavior  of  the  complete  rocket  engine  have  been  initiated  with  a  large 

amplitude  disturbance.  Since  longitudinal  combustion  instabilities  in 

actual  engines  seem  to  arise  spontaneously,  they  presumably  originate 

from  small  amplitude  unorganized  disturbances  which  may  be  present  in 

*  . 

the  chamber.  The  propellant  system  with  Q  =  -  120  cal/gm,  whose  com- 

s 

bustion  response  curve  is  shown  in  Figure  42,  was  employed  in  a  rocket 
engine  [see  Case  (9)5  Appendix  B]  designed  to  operate  near  the  peak 
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response.  The  propellant  mass  flow  rate  should  be  in-phase  with  a  small 
amplitude  pressure  oscillation  at  the  fundamental  frequency  of  the 
chamber  in  Case  (9).  The  steady-state  flow  field  for  this  motor  was 
altered  with  a  ten  percent  amplitude  continuous  disturbance,  with  the 
same  shape  as  shown  in  Figure  33*  The  time-history  of  the  pressure  at 
the  head  end  of  the  chamber  is  shown  in  Figure  48  for  the  restarted 
computation.  The  instantaneous  propellant  mass  flow  rate  at  the  same 
location  is  also  shown  in  the  figure.  The  results  show  that  the  dis¬ 
turbance  is  sustained  at  approximately  the  initial  amplitude  for  two 
cycles.  However,  after  the  second  reflection,  the  mass  flow  rate  spikes 
creating  a  related  jump  in  pressure.  As  mentioned  before,  the  trun¬ 
cation  error  introduced  into  the  time -dependent  flow  field  by  this 
behavior  of  the  propellant  burning  rate  will  compromise  its  accuracy. 
Thus,  the  flow  field  response  beyond  t  =  1.3  must  be  interpreted  as  a 
trend.  The  results  show  an  unmistakable  trend  toward  a  large  amplitude 
disturbance  which  becomes  a  shock  wave  by  the  sixth  reflection.  Thus, 
within  the  accuracy  limits  of  the  present  computational  scheme,  a  pro¬ 
pellant  system  described  by  typical  values  of  an  AP-propellant  is  shown 
to  be  capable  of  driving  a  small  amplitude  disturbance  into  a  large 
amplitude  axial  instability  in  the  solid  propellant  rocket  engine. 


CHAPTER  V 


CONCLUSIONS  AND  RECOMMENDATIONS 

A  mathematical  model  describing  the  behavior  of  a  solid  propellant 
rocket  engine  with  a  cylindrically  perforated  propellant  grain  has  been 
developed  to  predict  the  time-dependent  response  of  the  engine  flow 
field  to  arbitrary  disturbances.  The  solution  includes  the  unsteady 
flow  in  the  nozzle.  The  equations  describing  the  pressure-coupled 
combustion  response  of  the  burning  solid  propellant  are  solved  simul¬ 
taneously  with  the  chamber  and  nozzle  flow  field.  The  results  obtained 
confirm  that: 

(a)  The  unsteady  flow  in  the  choked  nozzle  is  important  in 
determining  how  disturbances  in  the  engine  are  reflected  in  the  throat 
region.  Specifically,  when  large  amplitude  disturbances  are  present, 
the  Mach  number  in  the  geometric  throat  cannot  be  assumed  to  remain  at 
unity. 

(b)  The  ability  of  typical  diameter,  inert,  mono-disperse  solid 
particles  to  damp  low  frequency  axial  motion  is  minimal.  Although  the 
drag  force  due  to  the  suspended  inert  particles  will  smooth  out  higher 
frequency  disturbances,  the  overall  effect  of  inert  particles  is  secon¬ 
dary.  However,  if  properly  modeled,  the  actual  particle  formation  pro¬ 
cess  and  the  energy  release  associated  with  metal  combustion  in  the  gas- 
phase  may  play  a  dominant  role  in  longitudinal  instability,  under  cer¬ 
tain  conditions . 
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(c)  For  all  flow  fields  examined,  the  steady  state  propellant 
burning  rate  law,  r  =  apn(n  <  l),  was  incapable  of  sustaining  longi¬ 
tudinal  instabilities.  This  indicates  that  the  dynamic  burning  rate 
characteristics  of  solid  propellants  cannot  be  ignored  in  a  study  of 
the  transient  behavior  of  a  solid  propellant  rocket  motor. 

The  following  conclusions  can  be  drawn  from  the  present  results 
for  rocket  engine  flow  field  response  to  arbitrary  disturbances: 

(1)  A  pressure-coupled  solid  propellant  combustion  model,  based 
on  commonly  employed  assumptions,  predicts  that  large  amplitude  burning 
rate  spikes  will  occur  under  certain  conditions  when  the  net  heat 
release  in  the  surface  reaction  is  exothermic.  Thus,  the  details  of 
the  condensed  phase  of  the  combustion  process  near  the  surface  are 
extremely  important . 

The  predictions  of  this  model  must  be  interpreted  within  the 
constraints  of  the  assumptions.  Certainly,  including  a  more  sophisti¬ 
cated  model  of  condensed  phase  reactions  and  treating  the  gas-phase 
flame  as  truly  time-dependent  may  alter  these  predictions  somewhat.  The 
important  point  is,  however  that  small  amplitude  pressure  oscillations 
in  the  flame  do  not  necessarily  result  in  small  amplitude  propellant 
burning  rates . 

(2)  With  this  combustion  model,  the  numerical  solution  for 

the  rocket  engine  flow  field  predicts  the  possibility  that  small  ampli¬ 
tude  disturbances  may  grow  into  large  amplitude  axial  instabilities. 
These  numerical  results  must  be  interpreted  as  a  trend  when  the  burning 
rate  exhibits  spikes,  as  a  result  of  the  unavoidable  truncation  errors 
in  the  time  integration  scheme.  However,  it  appears  that  even  without 
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the  contribution  of  velocity  coupling,  the  combustion  response  due  to 
pressure  coupling  alone  has  the  potential  for  driving  large  amplitude 
axial  instabilities. 

(3)  Marxman's  conclusion  that  a  combustion  response  magni¬ 

tude  of  3  to  4  is  the  minimum  value  needed  to  sustain  large  amplitude 
axial  instability  is  qualitatively  correct.  However,  the  phase  relation¬ 
ship  between  the  mass  flow  rate  of  the  burning  propellant  and  the  pres¬ 
sure  disturbance  is  extremely  important.  If  the  disturbances  in  the 
eigine  achieve  large  amplitudes,  a  propellarc  system  operating  on  the 
lead  portion  of  the  combustion  response  curve  will  be  considerably 
more  effective  in  driving  an  instability  than  one  operating  on  the  lag 
portion  of  the  curve. 

Toward  the  goal  of  making  quantitative  predictions  of  the  sus¬ 
ceptibility  of  a  given  rocket  engine  and  propellant  system  to  axial 
instability,  it  is  recommended  that  the  following  two  areas  receive 
further  attention: 

(1 )  The  solid  propellant  combustion  process  must  be  modeled  in 
greater  detail  by  removing  the  assumptions  of  a  quasi-steady  flame  and 
no  condensed  phase  reactions.  A  thorough  study  of  the  time -dependent 
behavior  of  this  model  should  indicate  what  changes  occur  in  the  spikes 
found  by  the  present  analysis.  If  the  sharp,  finite-amplitude  burning 
rate  spikes  remain,  then  a  new  method  of  integrating  the  conservation 
equations  in  the  chamber  will  be  required. 

(2)  To  deal  with  metal-loaded  composite  propellants,  the  solid 
particle  formation  process  including  the  combustion  of  molten  metal 
agglomerates  entrained  in  the  gas  flow  must  be  accounted  for  in  the 
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model.  This  has  the  potential  of  being  a  significant  energy  source  for 
driving  axial  instabilities. 

These  improvements,  added  to  the  capabilities  of  the  present  analysis, 
should  be  a  significant  step  toward  "a  priori"  predictions  of  axial 
instabilities  in  rocket  systems  in  the  design  stage. 


APPENDIX  A 


DERIVATION  OF  EQUATION  (3-30),  CHAPTER  III 

(54) 

Following  Sneddon's  '  '  work,  consider  a  linear  first-order 

partial  differential  equation  of  the  form, 


3z  Sz 

P(x,y,z)  ^  +  Q(x,y,z)  ^  =  R(x,y,z) 


(A.l) 


The  solution  to  Equation  (A.l)  is  an  arbitrary  function  of  the  solutions 
to  the  characteristic  equations 


dx  _  dy  _  dz 
P  =  Q  =  R 


(A.2) 


or 


§  -  »/«  (A.3) 

§  =  P/4  (A.l*) 


To  form  the  analogy  with  Equations  (3.28,  3-29)>  let 
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Preceding  page  blank 


T  -•  z 

v  -*  x 

I 

i 

y  ■*  y 

J 

i 

then, 

|  §  -  */«  (A-5) 

j 

» 

33?  -  (A-6) 

Then  if, 

A 

R(T,  V,  y)  -  v 

Q(T,  y,  y)  =  1 

P(T,  v,  y)  =  rv'  +  ^  (T  -  Tn_1) 
the  characteristic  equations  [Equations  (A. 5,  A. 6)]  become 
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I 


T'  =  v 


•1  A  A  y.  1 

v*  =  rv  +  77  (T  -  T  ) 

At 


which  are  identical  to  Equations  (3.28,  3 -29) •  The  solution  to  this 
system  must  then  yield  the  solution  to 


dT  ,  BT  r  ,  1  /A  An-1*., 

+  to  [rv  +  At  (T  "  T  »  "  T 


by  comparison  with  Equation  (A.l).  This  is  Equation  (3  - 30) . 
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APPENDIX  B 


ENGINE  AND  PROPELIANT  PARAMETERS  USED  IN  COMPUTATIONS 


Case  (l) 

Chamber  radius  =  1.0  ft. 
Engine  length  =  10.0  ft. 
Grain  length  =  8.0  ft. 
Nozzle  area  ratio  =4.0 


Chamber  temperature 
Chamber  pressure 
m 

g 

Isentropic  index 


4500°R 
1950  psia 
O.lp0-3 
1.20 


Case  (2) 


Chamber  temperature 
Chamber  pressure 


(a)  without  particles: 

Chamber  radius  =  1.0  ft. 

Engine  length  =  10.0  ft. 

Grain  length  =  8.0  ft. 

Nozzle  area  ratio  =  11.65  Isentropic  index 

(b)  vith  10fx  diameter  paricles 


in 


g 


Chamber  radius  =  1.0  ft. 
Engine  length  =  10.0  ft. 
Grain  length  =  8.0  ft. 


Chamber  temperature 
Chamber  pressure 


m 


g 


Nozzle  area  ratio  =11.65  m 


K 


-  54.56 


Isentropic  index 


4500°R 
397  psia 

0.008p0,1* 

1.20 

4500°R 
394  psia 

0.008p0'^ 

0.1  m 

g 

1.20 


Case  (3) 

Flame  temperature  =  4500° R 

Flame  pressure  =  1000  psia 

Isentropic  index  =  1.20 
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Preceding  page  blank 


Mean 

molecular  weight  = 

29  lb  /lb  -mole 
'  nr  m 

* 

Qs  = 

-100  cal/gra 

11 

545  cal/gm 

* 

c  = 
s 

0.275  cal/gm-°K 

* 

cp  = 

0.300  cal/gm-0 K 

* 

k  = 
s 

1.2  x  10”')  cal/cm- 

11 

1*  to 

_k 

2.0  x  10  cal/cm-sec-°K 

sec-°K 

* 

ps  = 

1-95  gm/cm3 

* 

Ef  = 

=  30  kcal/mole 

* 

E  = 
s 

20  kcal/mole 

*'r  *  <£ 

i-'  =  1.0  cm/sec  @  p  =  1000  psia,  Tsur  =  88o°K 

All  the  cases  listed  below  include  the  following  values  unless  other¬ 
wise  noted: 

Solid  propellant: 

p*  =1.95  gm/cm3  =  121.8  lbm/ft3 

T*  =  300°  K  =  54o°k 
cs 

c*  =  0.275  cal/gm-°K  a  0.275  BTU/lbm-°R 
k*  =  1.20  X  10"3  cal/cra-sec-°K  =  8.06  x  10~5  BTU/ft-sec-°R 
Eg  =  20  kcal/mole 
Gas  phase: 

V  =  1.20 

T*  =  2500°K  =  4500° R 

c*  =  0.300  cal/gm-°K  =  0.300  BTU/lbm-°R 

•X-  .ill 

k  =  2.0  x  10  cal/cm-sec-°K 
g  ' 

=  1.342  x  10"5  BTU/ft-sec-°R 
E^  =  30  kcal/mole 

W  =22  gm/gm-mole  =  22  lb  /lb  -mole 
u  nr  m 
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Engine: 


Chamber  radius  =  1.0  ft. 

Engine  length  =  10.0  ft. 

Grain  length  =  8.0  ft. 

Case  (4) 

Nozzle  area  ratio  =  11.65 
Chamber  pressure  =  396  psia 
Q*  =  1053.  BTU/lb^  (endothermic) 

Q*  =  2215.  BTU/lbm  (exothermic) 

r*  =  0.507  cm/sec  @  p  =  400  psia,  T  =  880°K 

'  SUT 


Case  (5) 

Nozzle  area  ratio  =  11. 65 
Chamber  pressure  =  397  psia 
Q*  =  -  180.0  BTU/lb^  (exothermic) 

Q*  =  981.8  BTU/lbm  (exothermic) 

r*  =  0.507  cm/sec  @  p*  =  400  psia,  Tgur  =  880°K 

Case  (6) 

Nozzle  area  ratio  =11.1 
Engine  length  =  15. 0  ft. 

Grain  length  =  13.0  ft. 

Chamber  pressure  =  1010  psia 
Q*  =  -  180.0  BTU/lb^  (exothermic) 

13b 


Q*  =  981.8  BTU/lbm  (exothermic) 

* 

r  =  1  cra/sec  @  p  =  1000  psia,  T  =  880°K 

Case  (7) 

Nozzle  area  ratio  =  11.65 
Chamber  pressure  =  397  psia 
Q*  =  -  207.0  BTU/lbm  (exothermic) 

Qf  =  95^.8  BTU/lb^  (exothermic) 

r  =  0.507  cm/sec  @  p*  =  400  psia,  T  =  880°K 
'  sur 


Case  (8) 

Nozzle  area  ratio  =  11.65 

Chamber  pressure  =  397  psia 

Q*  =  -  220.0  BTU/lb  (exothermic) 
s  m 

=  941.8  BTU/lb^  (exothermic) 

r  =  0.507  cm/sec  @  p  =  400  psia,  =  880°K 


Case  (9) 

Nozzle  area  ratio  =  12.1 
Grain  length  =  8.5  ft. 
Chamber  pressure  =  600  psia 

-g. 

Qg  =  -  216.0  BTU/lb^  (exothermic) 
0^  =  945.8  BTU/lb^  (exothermic) 


14u 


=  0.685  cm/sec 


600  psia, 


# 

T 

sur 


880°K 
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APPENDIX  C 


SAMPLE  CALCULATION  OF  CONSTANTS  IN  COMBUSTION  THEORY 

The  following  computation  applies  to  Case  (4)  of  Appendix  B 
where  the  heat  release  associated  with  the  surface  reaction  was  assumed 
to  be  1,053  BTU/lb  ,  endothermic.  The  information  below  is  assumed 
known  from  experimental  data: 

r  =  1.667  x  10  2  ft/sec 

at  p  =  400  psia,  T  =  1585°R 

*  *  ’  sur 

This  is  in  addition  to: 

p*  =  121.8  lb  /ft3 
rs  m' 

c*  =  0.275  BTU/lb  °R 
s  '  m 

k*  =  8.06  x  10"5  BTU/ft-sec-°R 

Eg  =  20  kcal/mole 

T*  =  54o°r 
cs 

c*  =  0.300  BTU/lb  °R 
p  '  m 

k*  =  1.342  x  10"5  BTU/ft-sec-°R 
T*  =  4500°R 
Eg  =  30  kcal/mole 

Y  =  1.20 


Preceding  page  blank 
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W  =  22  lb  /lb  -mole 
m'  m 

*  -X* 

Three  parameters  (B  ,  Q_,  m  )  must  be  determined. 
r  r  i  con 


*  *  Es/RoTsur 

r  =  B  e 
r 


Equation  (3*3) 


* .  * 

.  ,  +E  /R  T  - 

-*  B*  =  r  e  8  0  sur  =  1.535  x  lO-3  ft/sec 


(2)  Overall  heat  balance  yields, 


*  *  #.  *  #  .  *.  #  *  . 

Q  =  Q  +  c  T  -  T  )  +  c  (T_  -  T  ) 
f  s  s  sur  cs  p'  f  sur' 


=  2,215  BTU/lb 


(3)  Integrated  energy  equation  applied  across  interface, 


rk*^Ll 

s  Sy*  0 


*1  =  -  mV  +  [k*  ^1 

*1  n  s  L  g  **] 


By  J  © 


where 


Tyl  0  ^  “  Zlr  +  Z2/r  Equation  (3-19) 


*  .  *  *  *  ,  *A 

h  5  VcsTr  >  r  *  r  /Vc 


144 


Now,  at  steady  state 


T  I  ^  =  (T  -  T  )r 
y'  0  '  sur  cs' 


Rearranging  Equation  (3-19) 


r  +  Z^2) 
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=  3.159  x  3.°4  lbm/ft3-sec 

if  L*  =  10  ft,  ar*  =  1,208  ft/sec  T*  =  540°R 

*  .  „ 
p  =14.7  psia 


s 


s 

t 
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APPENDIX  D 

FIGURES  FOR  CHAP1ER  IV 
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pressure,  p  (atmospheres) 
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time,  t  (non-dimensional) 

Figure  15.  Pressure-Time  History  at  the  Head  End  of  the  Combustion 
Chamber  During  the  Ignition  Transient  For  Case  (l) 
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Figure  l6.  Pressure  Distributions  in  Rocket  Engine  Flow  Field 
of  Case  (l)  as  the  Result  of  a  Forty  Percent 
Amplitude  Shock  Wave  Disturbance 


pressure,  p  (atmospheres) 
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Figure  17.  Pressure-Time  History  at  the  Head  End  of  the 
Combustion  Chamber  as  the  Result  of  a  Forty 
Percent  Amplitude  Shock  Wave  Disturbance  on 
the  Steady  State  Flow  Field  of  Case  (l) 
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Figure  20.  Pressure-Time  History  at  the  Head  End  of  the  Combustion 
Chamber  as  the  Result  of  a  Forty  Percent  Amplitude 
Continuous  Disturbance  Added  to  the  Steady  State  Flow 
Field  of  Case  (2) 
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6.0 


Low  Frequency  (tu  =  2tt) 


Figure  22.  Propellant  Surface  Mass  Flow  Rate  Due 
to  a  Ten  Percent  Amplitude  Sinusoidal 
Pressure  Oscillation  at  a  Low  Frequency 
(u)  =  2rt)  and  a  High  Frequency  (u)  =  48tt) 
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Figure  2k.  Comparison  of  Combustion  F asponse  to  a  Ten  Percent  Amplitude  Square 
Wave  Pressure  Oscillation,  and  a  Sine  Wave  Pressure  Oscillation  at 
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figure  28.  Time  History  of  Pressure  and  Propellant  Mass  Flow  Rate  at  the 
Head  End  of  the  Combustion  Chamber  as  the  Result  of  a  Forty 
Percent  Amplitude  Shock  Wave  Disturbance  on  Case  (5) 


shock 


Figure  29.  Time  History  of  Pressure  and  Propellant  Mass  Flow  Rate  at  the 
Head  End  of  the  Combustion  Chamber  as  the  Result  of  a  Forty 
Percent  Amplitude  Shock  Wave  Disturbance  on  Case  (6) 
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Figure  30.  Combustion  Response  as  a  Function  of  Frequency 
fop.  the  Propellant  System  in  Case  (7) 
lQ„  =  -115  cal/gm] 
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Figure  32.  Combustion  Response  as  a  Function  of  Frequency  for  the 
Propellant  System  in  Case  (8)  [Q*  =  -  .  ?P  cal/gm] 


pressure,  p  (atmospheres) 


Figure  33.  Continuous  Type  Disturbances  Added  to  the  Steady  State 
Flow  Field  of  Case  (8) 
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Figure  3^-  Time  History  of  Pressure  and  Propellant  Mass  Flow  Rate  at  the 
Head  End  of  the  Combustion  Chamber  as  the  Result  of  a  Twenty 
Percent  Amplitude  Continuous  Disturbance  on  Case  (8) 


time,  t  (non-dimensional) 


Figure  35.  Time  History  of  Pressure  and  Propellant  Mass 
Flow  Rate  at  the  Head  End  of  the  Combustion 
Chamber  as  the  Result  of  a  Forty  Percent 
Amplitude  Continuous  Disturbance  on  Case  (8) 
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uonoustion  Kesponse  to  a  Five  Percent  Amplitude  Sinusoidal 
Pressure  Oscillation  at  0  =  10.0  for  the  Propellant  System 
in  Case  (8) 
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Figure  37.  Combustion  Response  to  a  Ten  Percent  Amplitude  Sinusoidal 
Pressure  Oscillation  at  f)  =  10.0  for  the  Propellant  System 
in  Case  (8) 


Figure  39-  Combustion  Response  to  a  Twenty  Percent  Amplitude  Sinusoidal 
Pressure  Oscillation  at  n  =  10.0  for  the  Propellant  System 
in  Case  (8) 
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Sinusoidol  Pressure  Oscillation  m  Flame  (1%  Shown) 


Figure  43.  Combustion  Response  to  Varying  Amplitude 
Sinusoidal  Pressure  Oscillations  at  0  = 
5>5  for  the  Propellant  System  in  Case  (9) 
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t  (non-dimensional) 


O.5I13  mil-sec 


Expanded  View  of  "Spike"  in  Pro£cl] ant 
Burning  Rate  Occurring  Between  t  -  3.2 
and  t  ~  3.30  in  Figure  43d 


4.00  (1200° K) 


Figure  45.  Temperature  Distributions  in  Unburned  Solid  Propellant  During 
Burning  Rate_Spike  Shown  in  Figure_44:  1,  Before  (t  =  3-088); 
2,  At  Peak  (t  =  3-245);  3,  After  (t  =  3-283) 
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Figure  b6.  Combustion  "espouse  to  a  Ten  Percent  Amplitude  Sinusoidal 
Pressure  Oscillation  at  n  =  11.0  for  the  Propellant  System 
in  Case  (9) 
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(9),  Altered  with  E*  =  10  kcal/mole 


pressure,  p  (atmospheres) 
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